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^ Abstract 

This survey concerns the study of quasi-stationary distributions with a specific 
Q focus on models derived from ecology and population dynamics. We are concerned 

Q with the long time behavior of different stochastic population size processes when 

is an absorbing point almost surely attained by the process. The hitting time of 
this point, namely the extinction time, can be large compared to the physical time 

and the population size can fluctuate for large amount of time before extinction 

actually occurs. This phenomenon can be understood by the study of quasi-limiting 
p I distributions. In this paper, general results on quasi-stationarity are given and 

(-H examples developed in detail. One shows in particular how this notion is related 

to the spectral properties of the semi-group of the process killed at 0. Then we 
g study different stochastic population models including nonlinear terms modeling 

i_i the regulation of the population. These models will take values in countable sets 

^ (as birth and death processes) or in continuous spaces (as logistic Feller diffusion 

processes or stochastic Lotka-Volterra processes). In all these situations we study 
CN in detail the quasi-stationarity properties. We also develop an algorithm based on 

Fleming- Viot particle systems and show a lot of numerical pictures. 

Keywords: population dynamics, quasi-stationarity, Yaglom limit, birth and death 
I process, logistic Feller diffusion, Fleming- Viot particle system. 
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1 Introduction 

We are interested in the long time behavior of isolated biological populations with a 
regulated (density-dependent) reproduction. Competition for limited resources impedes 
these natural populations without immigration to grow indefinitely and leads them to 
become extinct. When the population's size attains zero, nothing happens anymore and 
this population's size process stays at zero. This point is thus an absorbing point for the 
process. Nevertheless, the time of extinction can be large compared to the individual time 
scale and it is common that population sizes fluctuate for large amount of time before 
extinction actually occurs. For example, it has been observed in populations of endangered 
species, as the Arizona ridge-nose rattlesnakes studied in Renault-Ferriere-Porter [52j , that 
the statistics of some biological traits seem to stabilize. Another stabilization phenomenon 
is given by the mortality plateau. While demographers thought for a long time that 
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the rate of mortality of individuals grows as an exponential function of the age, it has 
been observed more recently that the rate of mortality slows at advanced ages, or even 
stabilizes. To capture these phenomena, we will study the long time behavior of the 
process conditioned on non extinction and the related notion of quasi-stationarity. In 
particular, we will see that a Markov process with extinction which possesses a quasi- 
stationary distribution has a mortality plateau. 

In all the following, the population's size process {Zt,t > 0) will be a Markov process 
going almost surely to extinction. We are interested in looking for characteristics of this 
process giving more information on its long time behavior. One way to approach this 
problem is to study the "quasi-limiting distribution" (QLD) of the process (if it exists), 
that is the limit, as t — )■ +00, of the distribution of Zt conditioned on non-absorption 
up to time t. This distribution, which is also called Yaglom limit, provides particularly 
useful information if the time scale of absorption is substantially larger than the one of 
the quasi-limiting distribution. In that case, the process relaxes to the quasi-limiting 
regime after a relatively short time, and then, after a much longer period, absorption will 
eventually occur. Thus the quasi-limiting distribution bridges the gap between the known 
behavior (extinction) and the unknown time-dependent behavior of the process. 
There is another point of view concerning quasi-stationarity. A quasi-stationary distribu- 
tion for the process [Zt, t > 0) denotes any proper initial distribution on the non-absorbing 
states such that the distribution of Z^ conditioned on non-extinction up to time t is inde- 
pendent of t, t > 0. 

There is a large literature on quasi-stationary distributions and Yaglom limits (see for 
example the large bibliography updated by PoUett [50]) and a lot of references will be 
given during the exposition. The present paper is by no means exhaustive, but is a sur- 
vey presenting a collection of tools for the study of QSD concerning specific population's 
size models. More than the originality of the proofs, we emphasize some general pat- 
terns for qualitative and quantitative results on QSD. We also provide a lot of numerical 
illustrations of the different notions. 

In Section [2] of this survey, we will introduce the different notions of QSD and review 
theoretical properties on QSD and QLD. We will also highlight the relations between 
QSDs and mortality plateaus. In Section [3| we will study the simple case of QSD for 
processes in continuous time with finite state space. We develop a simple example to 
make things more concrete. Thus we will concentrate on QSD for several stochastic 
population models corresponding to different space and time scales. We will underline 
the importance of spectral theory as mathematical tool for the research of QSD, in these 
different contexts. In Section |4[ we will consider birth and death processes. We will state 
results giving explicit conditions on the coefficients ensuring the almost sure extinction of 
the process, and the existence and uniqueness (or not) of a QSD. We will especially focus 
on the density-dependence case, when the death rate of each individual is proportional 
to the population's size (called logistic birth and death process). We will show that in 
that case, the process goes almost surely to extinction, and that there is a unique QSD, 
coinciding with the unique QLD. In Section [5| the birth and death process is rescaled 
by a growing initial population and by small individual weights (small biomass). This 
process is proved to converge, as the initial population's size tends to infinity, to the 
unique solution of the deterministic logistic equation, whose unique stable equilibrium is 
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given by the carrying capacity. If the individual birth and death rates are proportional 
to the population's size while preserving the ecological balance, more numerous are the 
individuals, smaller are their weights and faster are their birth and death events, reflecting 
allometric demographies. In that case, the rescaled birth and death process converges, as 
the initial population size increases, to the solution of a stochastic differential equation 
with a 1/2-Holder diffusion coefficient and a quadratic drift, called logistic Feller equation. 
The existence of the QSD is proved in this case and uniqueness is characterized by a 
condition meaning the return of the process from infinity in finite time. The proof relies 
QSD investigation to spectral theory tools developed in a functional L^— space. The 
logistic Feller equation describes the size of a mono-type population where individuals 
are indistinguable. Motivated by ecological and biodiversity problems, we generalize the 
model to multi-type populations with intra-specific and inter-specific competition. That 
leads us to consider stochastic Lotka-Volterra processes. We give conditions ensuring 
mono-type transient states or coexistence, preserving one or several dominant types in a 
longer scale. A large place in the survey is given to illustrations obtained by simulations. 
Some of them derive from an approximation method based on Fleming- Viot interacting 
particle systems, which is carefully described in Section [6j 

A brief bibliography on quasi-stationary distributions 

The study of quasi-stationary distributions began with the work of Yaglom on sub-critical 
Galton- Watson processes j6^. Since then, the existence, uniqueness and other properties 
of quasi- stationary distributions for various processes have been studied. 

In particular, the case of Markov processes on finite state spaces has been studied by 
Darroch and Seneta, who proved under some irreducibility conditions the existence and 
uniqueness of the QSD, for both discrete [TB] and continuous time settings [19j (detailed 
proofs and results are reproduced in Section [s] of the present paper). We also refer the 
reader to the works of van Doorn and Pollett |61| for a relaxation of the irreducibility 
condition. 

The case of discrete time birth and death processes has been treated by Seneta and 
Vere- Jones |54] and Ferrari, Martinez and Picco j2lj. For continuous time birth and 
death processes, we refer to van Doorn [59]. This last case is quite enlightening, since 
it leads to examples of processes with no QSD, of processes with a Yaglom limit and 
an infinite number of QSD and to processes with a Yaglom limit which is the unique 
QSD (detailed proofs and results are also developed in Section |4] of this survey). For 
further developments, we may refer to Pakes and Pollett |48| (where results on continuous- 
time birth and death processes with catastrophic events are obtained), to Bansaye [5] 
(where a discrete time branching process in random environment is studied), to Coolen- 
Schrijner [T7] (where general discrete time processes are studied) and references therein. 

Diffusion processes have also been extensively studied in the past decades, beginning 
with the seminal work of Mandl for the one-dimensional case and of Pinsky [33] and 
Gong, Qian and Zhao [29] in the mult i- dimensional situation. Martinez, Picco and San 
Martin |16] and Lladser and San Martin |44] highlighted cases of diffusions with infinitely 
many quasi-limiting distributions, with a non-trivial dependence on the initial distribution 
of the process. For recent development of the theory of QSDs for diffusion processes, we 
refer to Steinsaltz and Evans [57j and Kolb and Steinsaltz [4Tj where the case of one 
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dimensional diffusions with different boundary conditions is studied. We also emphasize 
that in the case of Wright-Fisher diffusions and some of its relatives, Huillet [35j derived 
explicit values of QSDs. Other diffusion processes related to demographic models have 
been studied in Cattiaux, Collet, Lambert, Martinez, Meleard and San Martin [13], where 
the case of the Feller logistic diffusion is developed (proofs and results are also written 
in detail in Section |5] of this paper), and in Cattiaux and Meleard [14j, where the case 
of a two dimensional stochastic Lotka-Volterra system is developed {k types stochastic 



Lotka-Volterra systems are also studied in Section 5.4 of this survey). 

Let us mention the original renewal approach of Ferrari, Kesten, Martinez and Picco [22j, 
also studied recently by Barbour and Pollett [6j in order to provide an approximation 
method for the QSD of discrete state space Markov processes in continuous time. Other 
approximation methods have been proposed by Pollett and Stewart [51] and by Hart and 
Pollett [34j . In this survey, we describe the approximation method based on Fleming- Viot 
type interacting particle systems, introduced by Burdzy, Holyst, Ingerman et March [lOj 
in 1996 and studied later by Burdzy, Holyst and March [TT], Grigorescu and Kang [32j . 
Ferrari and Marie [23], Villemonais [65] [66] and Asselah, Ferrari and Groisman [3]. 

For studies on the so-called Q-process, which is the process distributed as the original 
process conditioned to never extinct, we refer the reader to the above cited articles [45J, 
[49], [29], [IH], [in] and, for further developments, to the works of Collet, Martinez and 
San Martin [16] and of Lambert [43] and references therein. 

The framework 

Let us now introduce our framework in more details. The population's size {Zt,t > 0) 
is a Markov process taking values in a subset of N or M^., in a discrete or continuous 
time setting. If the population is isolated, namely without immigration, then the state 0, 
which describes the extinction of the population, is a trap. Indeed, if there are no more 
individuals, no reproduction can occur and the population disappears. Thus if the system 
reaches 0, it stays there forever, that is, if = for some t, then Zg = for any s >t. 

We denote by Tq the extinction time, i.e. the stopping time 

To = mi{t>0,Zt = 0}. (1) 

We will consider cases for which the process goes almost surely to zero, whatever the 
initial state is, namely, for all z G -E, 

P,(To < oo) = 1. (2) 

Before extinction, the process takes its values in the space E* = E\{0}. Any long time 
distribution of the process conditioned on non-extinction will be supported by E*. 

Notations For any probability measure n on E*, we denote by (resp. E^) the prob- 
ability (resp. the expectation) associated with the process Z initially distributed with 
respect to /i. For any x & E*, we set P^; = P^^ and E^^ = E^^. We denote by {Pt)t>o the 
semi-group of the process Z killed at 0. More precisely, for any z > and / measurable 
and bounded on E*, one defines 

PJiz)=E,if{Z,)h^To)- (3) 
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For any finite measure ^ and any bounded measurable function /, we set 

/^(/) = / fix)fx{dx), 

J E* 

and we also define the finite measure ^Pt by 

/xPi(/)=/i(^J)=E^(/(^t)l*<To). 

2 Definitions, general properties and first examples 

There are several natural questions associated with this situation. 

Question 1 What is the distribution of the size of a non-extinct population at a large 
time t ? The mathematical quantity of interest is thus the conditional distribution of Zt 
defined, for any Borel subset A d E* ^hy 

where v is the initial distribution of the population's size Zq. We will study the asymptotic 
behavior of this conditional probability when t tends to infinity. The first definition that 
we introduce concerns the existence of a limiting conditional distribution. 

Definition 1. Let a he a probability measure on E* . We say that a is a quasi-limiting 

distribution (QLD) for Z, if there exists a probability measure v on E* such that, for 
any measurable set A G E* , 

lim {Zt e A\To >t) = a{A). 

In some cases the long time behavior of the conditioned distribution can be proved to be 
initial state independent. This leads to the following definition. 

Definition 2. We say that Z has a Yaglom limit if there exists a probability measure 
a on E* such that, for any x E E* and any measurable set A G E* , 

\im¥^{ZteA\TQ>t) = a{A). (5) 

When it exists, the Yaglom limit is a QLD. The reverse isn't true in general and ^ will 
actually not imply the same property for any initial distribution. 

Question 2 As in the ergodic case, we can ask if this Yaglom limit has the conditional 
stationarity property given by the following definition. 

Definition 3. Let a be a probability measure on E* . We say that a is a quasi- stationary 

distribution (QSD) if, for all t > and any measurable set A G E* , 

a{A) = P„ (Zt e A\To > t) . 
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The main questions are: Does a QSD exists? Is there a unique QSD for the process? 
We will study examples where QSDs do not exist, or with an infinity of QSDs, or with 
a unique QSD. The relation between the existence of QSD, QLD and Yaglom limit is 
clarified in Proposition [T] below. Namely, we will prove that 

Yaglom limit QSD <^ QLD. 



Question 3 Since the processes we are interested in become extinct in finite time almost 
surely, the event t < Tq becomes a rare event when t becomes large. An important 
question is then to know whether the convergence to the Yaglom limit happens before the 
typical time of extinction, or if it happens only after very large time periods, in which 
case the populations whose size are distributed with respect to the Yaglom limit are very 



rare. Both situations can appear, as illustrated by the simple example of Section 2.3 

Question 4 While most of theoretical results on QLDs, QSDs and Yaglom limits are 
concerned with existence and uniqueness problems, it would be useful in practice to have 
qualitative information on the Yaglom limit. We present here particle approximation 
results and numerical computations of the Yaglom limit for some population's size models, 
providing some enlightenment on Question 3 above. 

Question 5 Another mathematical quantity related to this conditioning is based on a 
pathwise point of view. In the finite state space case of Section |3] and the logistic Feller 
diffusion case of Section |5j we will describe the distribution of the trajectories who never 
attain the trap. This will allow us to define a process, commonly referred to as the Q 
process for Z. We will prove that the new process defined by this distribution is ergodic, 
and that its stationary distribution is absolutely continuous with respect to the QSD (but 
not equal). 



The present section is organized as follows. In Subsection 2.1, we state general properties 



of QLDs, QSDs and Yaglom limits. In Subsection 2.2, we develop the case of the Galton- 
Watson process. This discrete time process is of historical importance, since the notion 
of Yaglom limit has originally been developed for this process by Yaglom itself (see [67j). 



In Subsection 2^, we develop a very simple example of a process evolving in a finite 
subset of N. For this process, one can easily prove the existence of the Yaglom limit, the 
uniqueness of the QSD, and compare the speed of extinction to the speed of convergence 
to the Yaglom limit. We also provide numerical computation of the relevant quantities. 



2.1 General properties 

Most of the following results are already known by the QSD community. In this section, 
we emphasize their generality. 



2.1.1 QSD, QLD and Yaglom limit 



It is clear that any Yaglom limit and any QSD is also a QLD. The reverse implication has 
been proved by Vere- Jones |63] for continuous time Markov chains evolving in a countable 
state space. The following proposition extends this result to the general setting. 
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Proposition 1. Let a be a probability measure on E* . The distribution a is a QLD for 
Z if and only if it is a QSD for Z . 

Remark 1. When it exists, tlie Yaglom limit is uniquely defined, while there are pro- 
cesses with an infinity of QSDs (see the birth and death process case of Section |4]). We 
immediately deduce that there exist QSDs which aren't a Yaglom limit. 

Proof. (1) If a is a QSD then it is a QLD for Z starting with distribution a. 

(2) Assume now that a is a QLD for Z and for an initial probability measure fi on E*. 
Thus, for any measurable and bounded function f on E*, 



a{f) = limE^(/(ZO|To>t)= lim 



E^(/(Zi);To>t) 



Applying the latter with f{z) = P2(To > s), we get by the Markov property 



Pc.(To > s) = lim 



P^(To >t + s) 



Let us now consider f{z) = F^^Zg G A, Tq > s), with A C E*. Applying the Markov 
property again, it yields 



Fa{Zs e A]To > s) = lim 



F^{Zt+seA;To>t + s) 



t^^ P^(To > t) 

^.^ Ff,{Zt+s eA;To>t + s) P^(To > t + s) 
P^(To>t + s) P^(To>t) 

By definition of the QLD a, ^^p^l'Ti^fl^)'^^^ converges to a{A) and ^'^^'^^j^^^f^^ converges 
to Pq,(To > s), when t goes to infinity. We deduce that, for any Borel set A of E* and 
any s > 0, 

aiA) = P«(Z, G A\To > s). 
The probability measure a is then a QSD. □ 

2.1.2 Exponential extinction rate 

Proposition 2. Let us consider a Markov process Z with absorbing point satisfying 
Assume that a is a QSD for the process. Then there exists a positive real number 9{a) 
depending on the QSD such that 

PjTo > t) = e-'^(")*. (6) 

This theorem shows us that starting from a QSD, the extinction time has an exponential 
distribution with parameter 6 {a) independent of t > 0, given by 

lnP^(ro>t) 
U[a) = . 
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Proof. By the Markov property, 

= P„(To > t)E, (Pz,(To > s)\To > t) , 

since Tq < t implies Zt = 0, and Po(To > s) = 0. By definition of a QSD, we get 

E„ {Fz,{To > s)\To >t)= P„(To > s). 

Hence we obtain tliat for all s,t > 0, Pa(To > t + s) = Pa(To > s)Pa(To > t). Let us 
denote g(t) = Po(To > t). We have 5f(0) = 1 and, because of ([2]), g{t) tends to as t tends 
to infinity. An elementary proof allows us to conclude that there exists a real number 
6 {a) > such that 

p^(To >t) = e-^(")*. 

□ 

2.1.3 QSD and exponential moments 

The following statement gives a necessary condition for the existence of QSDs in terms of 
existence of exponential moments of the hitting time Tq. 

Proposition 3. Assume that a is a QSD. Then, for any < 7 < 6{a), 

E«(e^^«) < +00. (7) 
In particular, there exists a positive number z such that E,z{e'^'^°) < +00. 

Proposition |3] suggests that if the population can escape extinction for too long times 
with positive probability, then the process has no QSD. This is the case for the critical 
Galton- Watson process: its extinction time is finite almost surely, but its expectation 
isn't finite. 

Proof. We compute the exponential moment in continuous and discrete time settings. In 
both cases, it is finite if and only if 6{a) > 7. 

In the continuous time setting, ^ says that, under P^,, Tq has an exponential distribution 
with parameter 6{a). We deduce that, for any 6{a) > 7, 

^ ' 0[a) — 7 

In the discrete time setting, ^ says that under P^, Tq has a geometric distribution with 
parameter e~^^"\ We deduce that 

Since E,a{e'^'^°) is equal to J^, Kz{e'^'^°)a{dz), the finiteness of the integral implies the last 
assertion. 

□ 
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Remark 2. In the particular case of an irreducible continuous time Markov chain with 
state space N such that lim2_j.+oo Pz (To < t) = 0, Vt > 0, Ferrari, Kesten, Martinez and 
Picco [22j proved that the existence of the moment ([T]) for some z G N and some 7 > is 
equivalent to the the existence of a quasi-stationary distribution. 

It is actually not true in any case, as shown by the following counter-example. Let Z 
be a continuous time random walk on N reflected on 1 and killed at rate 1. Thus, for 
any A G [0,1 [ and any probability measure /i on N, E^(e^'^°) is finite. Nevertheless the 
conditional distribution Fz{Zt G -jt < Tq) is the distribution of a standard continuous 
time random walk reflected on 1, which converges to as t tends to infinity. In particular 
Z has no QLD and thus no QSD. 



2.1.4 A spectral point of view 

In this section, the results are stated in the continuous time setting. The operator L with 
domain T'(L) denotes the infinitesimal generator of the sub-Markovian semi-group (Pt) 
associated with the killed process Z. The next proposition links the existence of QSDs 
for Z and the spectral properties of the dual of the operator L. It is one of the main tools 
used in a large literature studying QSDs. 

Proposition 4. Let a be a probability measure on E* . We assume that there exists a set 
D C "DlL) such that, for any measurable subset A G E* , there exists a uniformly bounded 
sequence (fn) in D converging point-wisely to 1^. 

Then a is a quasi- stationary distribution if and only if there exists 9 (a) > such that 

a{Lf) = -9{a)a{f), V/ G D. 

We emphasize that the existence of D is always true if the state space E* is discrete. 
It is also fulfilled if E* is an open subset of M"^ and if Z is a diffusion with locally bounded 
coefficients. 

Proof. (1) Let a be a QSD for Z. By definition of a QSD, we have, for every Borel set 
A C E*, 

a{A) = ; r. 



By Theorem |6| there exists 9{a) > such that for each t > 0, 

aPt{lE')=F^{To>t)=e-"^''^'. 



We deduce that, for any measurable set AGE*, aPt{lA) = e 
lent to aPf = e 
have 



a{A), which is equiva- 



By Kolmogorov's forward equation and by assumption on D, we 
dPtf, 



dt 



\PtLf{x)\ < IIL/IU <+oo, yfeD. 



In particular, one can differentiate aPtf = J^, Ptf{x)a{dx) under the integral sign, which 
implies that 

a{Lf) = ~9{a)aU), V/ G D. 



9 



(2) Assume know that a{Lf) = —9{a)a{f) for all f E D. By Kolmogorov's backward 
equation and the same "derivation under the integral sign" argument, we have 

= a{Lm = -9{a)aPt{f), V/ G D. 

We deduce that 

aPtif) = e-*^(")a(/), V/ G D. 

By assumption, there exists, for any measurable subset A G E*, & uniformly bounded 
sequence (/„) in D which converges point-wisely to 1^- Finally, we deduce by dominated 
convergence that 

This implies immediately that a is a quasi-stationary distribution for Z. □ 
2.1.5 Long time limit of the extinction rate 

Another quantity of interest in the demography and population's dynamics is given by 
the long time behavior of the killing or extinction rate. In the demography setting, the 
process Z models the vitality of some individual and t its physical age. Thus Tq is the 
death time of this individual. The long time behavior of the extinction rate has been 
studied in detail by Steinsaltz-Evans |56| for specific cases. 

The definition of the extinction rate depends on the time setting: 

• In the discrete time setting, the extinction rate of Z starting from /i at time t > 
is defined by 

=P^(To = t + l|To >t). 

• In the continuous time setting, the extinction rate of Z starting from fi at time t > 
is defined by 

' p,(ro > t) ' 

when the derivative exists and is integrable with respect to fi. 

Historically (c/. [28]). demographers applied the Gompertz law meaning that this extinc- 
tion rate was exponentially increasing with time. However in 1932, Greenwood and Irwin 
[31j observed that in some cases, this behavior was not true. In particular there exist 
cases where the extinction rate converges to a constant when time increases, leading to 
the notion of mortality plateau. This behavior of the extinction rate has been observed 
in experimental situations (see for instance [12]). 

The QSDs play a main role in this framework. Indeed, by Proposition |2| if a is a QSD, 
then the extinction rate r^it) is constant and given by 

_ f 1 — e~^^°^ in the discrete time setting ^ g 
" 1 9 (a) in the continuous time setting ' ~ 



We refer to the introduction of Steinsaltz-Evans [56] for a nice discussion of the notion of 
QSD in relationship with mortality plateaus. 

In the next proposition, we prove that the existence of a QLD for Z started from /i implies 
the existence of a long term mortality plateau. 
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Proposition 5. Let a be a QLD for Z , initially distributed with respect to a probability 
measure /i on E* . In the continuous time setting, we assume moreover that there exists 
h > such that L{PhlE*) is well defined and bounded. In both time settings, the rate of 
extinction converges in the long term: 

lim r^{t) = r„(0). (8) 

Proof. In the discrete time setting, by tlie semi-group property and tlie definition of a 
QLD, we liave 

r,{t) = l-^^^^^^ T— ^ l-a{P,lj,,) = r^{0). 

Tlie limit is by definition the extinction rate at time of Z starting from a. 

In the continuous time setting, by the Kolmogorov's backward equation, we have 

^^Pt+hlE^x) = PtL{PhlE*){x), Va; G E*. 

Since L{PhlE*) is assumed to be bounded, we deduce that 

d 



^^flPt+hi^E*) = jJ^PtHPhlE*)- 



Then 



§-fiPt+h{lE*) flPMPhlE*) _^ ,.p, ^ WP1 

- > a{LPhlE*) = -0{a)a{PhlE* 



by the definition of a QLD and by Proposition |4} We also have 

— — ^ ^ a{PhlE*)- 

Finally, we get 

r (t + h) = -^^^ — - — ^ > 0{a), 

fl{Pt+hlE*) 

which allows us to conclude the proof of Proposition |5} □ 



2.2 A historical example in discrete time: the Galton- Watson 
process 

The Galton- Watson process is a population's dynamics model in discrete time, whose size 
{Zn)n>o evolves according to the recurrence formula Zq and 

i=l 

where (^j-"'*)j,n is a family of independent random variables, identically distributed follow- 
ing the probability measure /i on N with generating function g. As defined, Zn is the 
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size of the n*^ generation of a population where each individual has a random number of 
children, chosen following yU and independently of the rest of the population. This process 
has been introduced by Galton and Watson (see [26j) in order to study the extinction of 
aristocratic surnames. 

We will assume in the whole section that < /i({0}) + < 1. We denote by 

m = E{^f'') the average number of children by individual in our Galton- Watson process. 
The independence of descendants implies that starting from Zq, the process Z is equal 
to the sum of Zq independent Galton- Watson processes issued from a single individual. 
By this branching property, the probability of extinction for the population starting from 
one individual is obtained as follows: 

Pi(3n G N, Z„ = 0) = lim Ei(0^") = lim ^ o ■ ■ ■ o ^(0) (n times). 

n—¥+oo n—>oo 

There are three different situations (see for instance Athreya-Ney 

- The sub-critical case m < 1: the process becomes extinct in finite time almost surely 
and the average extinction time E(To) is finite. 

- The critical case m = 1: the process becomes extinct in finite time almost surely, 
but E(To) = +00. 

- The super-critical case m > 1: the process is never extinct with a positive proba- 
bility, and it yields immediately that E(To) = +oo. 

Theorem 6 (Yaglom |67) . 1947). Let (Z„)„>o be a Galton-Watson process with the re- 
production generating function g. There is no quasi- stationary distribution in the critical 
and the super- critical case. In the sub-critical case, the Yaglom limit exists and is the 
unique QSD of Z . Moreover, its generating function g fulfills 

g{g{s)) = mg{s) + 1 - m, Vs G [0,1]. (9) 

Proof. The proof is adapted from Athreya-Ney [4j p. 13-14. In the critical or the super- 
critical case, we have Ei(To) = +oo, which implies that Ea(To) = +oo for all probability 
measure a on N*. We deduce from Proposition [3] that there is no QSD. 

Assume now that m < 1. Let us fix an arbitrary probability measure v on N* and 
prove that there exists a QLD a for Z starting with distribution v. For each n > 0, we 
denote by (?„ the generating function of Z„, gn{s) = E,^ (s^") Vs G [0,1]. Recall that 
gn+i = gi ° gn- Let us also denote by cjn the generating function of Z„ conditioned to 
{Zn > 0} = {To > n}: 



gn{s) = E,(s^"|Z„>0) 



P.(^n > 0) 

E.(g^")-P.(Z„ = 0) 
1 - P,(Z, = 0) 

1-gM l-9n{0) ^ ' ^' 
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Note that ^n(O) = 0, which is quite natural since the conditional law doesn't charge 0. 
For a fixed s G [0,1), we set T{s) = ^—^ij^- Then we have, for all n > 0, 

1 - 9n+l{s) = , , 1 - gn{s)) . 

r(^n(o)) 

Since gi is convex, F is non-decreasing. Moreover m < 1 implies that gn{x) > x, so that 
gn{s) and 1 — ^„(s) are non-decreasing in n. In particular, liuin^oo gn{s) exists. Let us 
denote by g{s) its limit and by a the corresponding finite measure (whose mass is smaller 
than one) . In order to prove that o; is a probability measure on N* , it is sufficient to prove 
that g{s) — i- 1 when s goes to 1. We have 

Tig^m (1 - ^„+i(s)) = (1 - gnigM)) . 

Taking the limit on each size, where lim„^oo r(5'n(0)) = F(l) = m, we deduce that 

m{l-g{s)) = 1 -g{gi{s)), 

which implies Equation Since lim<j_j.i (yfi(s) = 1 and m < 1, then ^(1) = 1. Finally, a 
is a QLD for Z starting with distribution u. 

One could think a priori that the function g depends on the starting distribution u. We 
prove now that it isn't the case, so that there is a unique QLD, and then a unique QSD, 
which is also the Yaglom limit of the process (indeed, one could choose u = 6x, x & W). 
Assume that there exist two generating functions g and h which fulfill Equation (|9]). By 
induction, we have, for all n > 1 and all s G [0,1], 

g{gn{s)) = m"^(s) + (m"-^ + ■ ■ ■ + m + l) (m - 1), 
h{gn{s)) = nf'hi^s) + (m""^ H h m + l) (m - 1). 

We deduce that for s G [0,1[ 

9\9n{s)) g'^s) = m-g'is) ; h'{gr,{s)) g'^{s) = m^h'^s). 

Since for the sub-critical case (?„(0) f 1 when n — )■ oo, for any s G [0,1[ there will be a A; 
such that 

9k{^)<s<gk+m. 

Hence, 

9'{s) _ g'jgnjs)) ^ g'{gn+k+im _ g'jO) mg'^^.jO) _ g'jO) m 
h'{s) h'{gn{s))~ h'{gn+km h'{<d) g'n+k+M h'{id) g'{9n+km' 

When n goes to infinity, we obtain < fyj^- The converse inequality is established 

similarly. Since g and h are generating functions of probability measures on N*, we have 
^(0) = /i(0) = and ^(1) = h{l) = 1. Finally, the two functions g and h are equal, which 
concludes the proof of Theorem |6] . □ 



13 



2.3 The simple example of an ergo die process with uniform killing 
in a finite state space 

Wc present a very simple Markov process with extinction whose quasi-stationary distri- 
bution, Yaglom hmit, speed of extinction and speed of convergence to the Yaglom hmit 

are very easy to obtain. 

Let {Xt) i>o be an exponentially ergodic Markov process which evolves in the state space 
E* — {1, ■ ■ ■ ,N}, N > 1. By exponentially ergodic, we mean that there exist a probability 
measure a on E* and two positive constants C,A > such that, for all z E {!,■ • ■ ,N} 
and all t >0, 

sup \F,{Xt = i)- a{{i})\ < Ce-^'. 
ieE* 

There is no possible extinction for (Xt). Let d > be a positive constant and let be 
an exponential random time of parameter d independent of the process (Xf). We define 
the process {Zt) by setting 

f Xt,ift<Td 

* \0, ift>r,. 

This model is a model for the size of a population which cannot be extinct, except at a 
catastrophic event which happens with rate d. Thus we have 

F^{t < To) = e""^*, Wt > 0. 

The conditional distribution of Zt is simply given by the distribution of Xf 

F,{Zt = i\Zt 7^ 0) = F,{Xt = i), Mz e E\ 

We deduce that the unique QSD is the Yaglom hmit a and that for all z e E* and all 
t > 0, 

sup \F^{Zt ^i\To>t)- a{{i})\ < Ce'^K 
ieE* 

Thus in this case, the conditional distribution of Z converges exponentially fast to the 
Yaglom limit a, with rate A > and the process becomes extinct exponentially fast, with 
rate d > 0. 

Hence the comparison between the speed of convergence to the Yaglom limit and the 
speed of extinction will impact the observables of the process before extinction: 

(a) If A ^ (i, then the convergence to the Yaglom limit happens before the typical time 
of extinction of the population and the quasi-stationary regime will be observable. 

(b) If A ^ (i, then the extinction of the population occurs before the quasi-stationary 
regime is reached. As a consequence, we are very unlikely to observe the Yaglom 
limit. 

(c) If A ~ (i, the answer is not so immediate and depends on other parameters, as in 
particular the initial distribution. 
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Survival Probability 




Distance to tlie Yaglom limit in the extinction's time scah 



Figure 1: Example 1. A numerical computation leads to A = 0.098. Three different situations are 
observed, which lead to three very different patterns for the speed of convergence to the Yaglom limit in 
the extinction's time scale: (o) A > d = 0.001; (□) A < d = 0.500; {■) X = d = 0.098. 



Example 1. The population size Z is described by a random walk in continuous time 
evolving in E = {0,1,2, ■ ■ ■ ,N} with transition rates given by 

z — )■ z + 1 with rate 1, for all i G {1,2, ■ ■ ■ — 1}, 
z — 7- z — 1 with rate 1, for all z G {2,3, ■ ■ ■ ,A^}, 
i with rate d> 0, for all i G {1,2, ■ ■ ■ ,A^}. 

The boundedness of the population size models a constraint of fixed resources which 
acts on the growth of the population. We will see more realistic fixed resources models 
including logistic death rate in the next sections. One can check that the quasi-stationary 
probability measure of Z is given hj ai = 1/N for all i & E*. 

Numerical simulations. We fix = 100. In that finite case, one can obtain by numerical 
computation the whole set of eigenvalues and eigenvectors of the infinitesimal generator 
L (we use here the software SCILAB). Numerical computation using the fact that A is 
the spectral gap of the generator of X gives A = 0.098. For different values d = 0.001, 
d = 0.500 and d = 0.098, we compute numerically the mathematical quantities of interest: 
the extinction probability P2(To > t) = e"'^* as a function of t (cf. Figure [l] left picture) 
and the distance sup^g^;. \Fz{Zt = z|To > t) — Q;({i})| between the conditional distribution 
of Zt and a as a function of — logF z{Tq > t), which gives the extinction's time scale, (cf. 
Figure [T] right picture). 



We observe that the convergence to the Yaglom limit happens rapidly in the case (o) A = 
0.098 ^ d = 0.001. Indeed the distance to the Yaglom limit is equal to 0.05, while 
the survival probability can't be graphically distinguished from 1. On the contrary, we 
observe that the convergence happens very slowly in the case (□) A = 0.098 -C d = 0.500. 
Indeed, the distance to the Yaglom limit is equal to 0.05 when the survival probability 
appears to be smaller than e^^^ ^ 3 x 10^''. The case (■) A = 0.98 = (i is an intermediate 
case, where the distance to the Yaglom limit is equal to 0.05 when the survival probability 
appears to be equal to ~ 0.05. 
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3 The finite case, with general kiUing's rate 



3.1 The quasi-stationary distributions 



The Markov process {Zt)t>o evolves in continuous time in E = {0,1,.. .,A^}, > 1 and 
we still assume that is its unique absorbing state. The semi-group {Pt)t>o is the sub- 
Markovian semi-group of the killed process and we still denote by L the associated in- 
finitesimal generator. In this finite state space case, the operators L and Pt are matrices, 
and a probability measure on the finite space E* is a vector of non-negative entries whose 
sum is equal to 1. The results of this section have been originally proved by Darroch and 
Seneta ([H] and [H]). 

Theorem 7. Assume that Z is an irreducible and aperiodic process before extinction, 
which means that there exists to > such that the matrix Pt^ has only positive entries (in 
particular, it implies that Pt has positive entries for t > to). Then the Yaglom limit a 
exists and is the unique QSD of the process Zt. 

Moreover, denoting by 6{a) the extinction rate associated to a (see Proposition^, 
there exists a probability measure n on E* such that, for any i,j G E* , 



lim e''("^*P,(Zt = j) 



TTi a 



t^oo 

and 

The main tool of the proof of Theorem [7] is the Perron-Frobenius Theorem, which 
gives us a complete description of the spectral properties of Pt and L. The main point is 
that the matrix Pi has positive entries. For the proof of the Perron-Frobenius Theorem, 
we refer to Gantmacher |27j or Serre [55j . 

Theorem 8 (Perron-Frobenius Theorem). Let (Pt) be a submarkovian semi-group on 
{1, • • • ,N} such that the entries of Pt^ are positive for to > 0. Thus, there exists a unique 
positive eigenvalue p, which is the maximum of the modulus of the eigenvalues, and there 
exists a unique left- eigenvector a such that > and Xlili ~ ^' '^^^ there exists a 
unique right-eigenvector vr such that VTj > and X^ili '^i'""* ~ ^' satisfying 

aPto = pa ] PtoTT = piT. (10) 

In addition, since (Pt) is a sub-Markovian semi-group, p < I and there exists 6 > such 
that p = . Therefore 

Pt = e-^'A + ^{e-^'), (11) 

where A is the matrix defined by Aij = TCiOj, and x > ^ '^i^'d "(9(6^^*) denotes a matrix 
such that none of the entries exceeds Ce~^*, for some constant C > 0. 

Proof of Theorem Applying Perron-Frobenius Theorem to the submarkovian semi-group 



(Pt)t>0; it is immediate from (11) that there exists 6 > and a probability measure a on 
E* such that, for any i,j G E*, 

e'%{Zt = j) = e''[Ptl, = -K.a, + ^(e-(^-^)*)- (12) 
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Summing over j E E* , we deduce that 

e^*P,(To>t)=7r, + ^9(e-('^-^)*). (13) 
It follows that, for any z,j G E*, 

Fi{Zt =j\To>t)= -— ^ aj. 

Fi{lo > t) t^oo 

Thus the Yaglom limit exists and is equal to a. Since E is finite, we have for any initial 
distribution u on E* 



lim F,{Zt = j|To > t) = V z/, hm Fi{Zt = j\To > t) = V] 
We deduce that the Yaglom limit a is the unique QLD of Z, and thus it is its unique 



QSD. By Proposition [2| we have aPi{lE*) = e By (10), this quantity is also equal 

to e^^, so that 9 = 9(a). The end of Theorem Mis thus a straightforward consequence of 



(12) and (13) □ 



Remark 3. One can deduce from (12) and (13) that there exists a positive constant Cl 
such that 

sup \Pi{Zt = j\Zt > 0) - a,| < CLe-(^-^(")), 
jeE*,ieE* 

where the quantity x ~ 9{a) is the spectral gap of L, i.e. the distance between the first 
and the second eigenvalue of L. Thus if the time-scale x ~ 9{a) of the convergence to 
the quasi-limiting distribution is substantially bigger than the time scale of absorption 
ix " 9{a) ^ 9(a)), the process will relax to the QSD after a relatively short time, and 
after a much longer period, extinction will occur. On the contrary, if x ~ 9{a) ^ 9{a), 
then the extinction happens before the process had time to relax to the quasi-limiting 
distribution. 

In intermediate cases, where A — 9{a) ~ 9{a), the constant Cl, which depends on the 
whole set of eigenvalues and eigenfunctions of L, plays a main role which needs further 
investigations. 



We generalize Example 1 to a more realistic case where the killing's rate can depend on 
the size of the population. For instance, it can be higher for a small population than for 
a big one. 

Example 2. Let Z he a Markov process which models a population whose individuals 
reproduce and die independently, with individual birth rate A > and individual death 
rate /i = 1. In order to take into account the finiteness of the resources, the process is 
reflected when it attains a given value N, that we choose here arbitrarily equal to 100. 
Thus the process Z evolves in the finite state space {0,1, ■ ■ ■ ,100} and its transition rates 
are given by 

z — 7- z + 1 with rate Xi, for all i G {1,2, ■ ■ ■ ,99}, 
z — 7- z — 1 with rate fii, for all z G {1,2,3, ■ ■ ■ ,100}. 
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(a) (b) (c) 



Figure 2: Example 2. Yaglom limits for different values of A. The following values of 9{a) are obtained 
by numerical computation, (a) A = 0.9, e{a) = 0.100; (b) A = 1.0, 0{a) = 0.014; (c) A = 1.1, 9{a) = 
5.84 X 10-5. 

The infinitesimal generator of Z is given by 
Li i = —1 — X and L12 = A, 

L,^i-i = i, Li^i = -{1 + \)i and Li^i+i = \i, Vi G {2, ■ ■ ■ ,99}, 

^100,99 = 100 and Lioo.ioo = -100, 

Lij = 0, Vi,j G {1, ■ ■ ■ ,100} such that \j - z| > 1. 

The process Z clearly fulfills the conditions of Theorem [7j As a consequence, it has a 
Yaglom limit a, which is its unique QSD. Moreover, the probability measure a is the 
unique normalized and positive left eigenvector of L. Since L is a finite matrix of size 
100 X 100, one can numerically compute the whole set of eigenvectors and eigenvalues of 
the matrix (Lij). This will allow to obtain numerically the Yaglom limit a, its associated 
extinction rate 6{a), and the speed of convergence x — 0{<y)- Moreover, for any t > 0, one 
can compute the value of e^^, which is equal to Pj (the semi-group of Z at time t). Hence, 
we may obtain the numerical value of the conditioned distribution FzdZt G .|t < Tq), 
for any initial size Zq. Finally, we are also able to compute numerically the distance 
between a and the conditioned distribution Fzo{Zt G .|t < Tq), for any value of A > and 
ZoE {!,■■■ ,100}. 

In Figure [2| we represent the Yaglom limit a for different values of A, namely A = 0.9, 
A = 1.0 and A = 1.1. Let us comment the numerical results. 

(a) In the first case (A = 0.9), an individual is more likely to die than to reproduce 
and we observe that the Yaglom limit is concentrated near the absorbing point 0. 
The rate of extinction 6{a) is the highest in this case, equal to 0.100. In fact, the 
process reaches the upper bound 100 very rarely, so that the behavior of the process 
is very similar to the one of a linear birth and death process with birth and death 
rates equal to A and fi respectively. In Section |4| we study such linear birth and 
death processes. We show that the Yaglom limit (which exists if and only if X < fi) 
is given by a geometric law and 6 (a) = fi — X. 

(b) In the second case (A = yU = 1), we observe that a decreases almost linearly from ai 
to aioo and the upper bound = 100 plays a crucial role. In fact, letting tend 
to +00, one would observe that for any i > 1, Oi decreases to 0. The extinction rate 
6 (a) which is equal to 0.014 for A^ = 100 would also go to 0. The counterpart of 
this phenomenon for the linear birth and death process studied in Section |4] is that 
the Yaglom limit will not exist when /i = A. 
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(c) Distance to the yaglom limit (c) Distance to the Yaglom limit in the extinction's time scale 

Figure 3: Example 2. Pictures (a) and (c) correspond to different values of A (the following values of 
^(q^) ~ X have been obtained by numerical computation): (a) A = 0.9, 0{a) — 0.100, 0{a) ~ x — 0.102; 
(c) A = 1.1, 9{a) = 5.84 x 10~^, 0{a) — x — 0.103; each curve corresponds to a given initial size of the 
population: (•) Zq = 1; (o) Zq = 10; (□) Zq = 100. 

(c) In the third case (A = 1.1), the Yaglom hmit a is concentrated near the upper 
bound 100, while the extinction rate is 6{a) = 5.84 x 10~^. The comparison with 
the linear birth and death process is no more relevant, since the important factor in 
this case is the effect of the upper bound = 100, which models the finiteness of 
the resources in the environment. 

In Figure [3} we study the effect of the initial position and of the value of the parameter 
A on the speed of convergence to the Yaglom limit and on the speed of extinction. We 
choose the positions Zq = 1, = 10 and Zq = 100, and we look at the two different cases 
A = 0.9 and A = 1.1, which correspond to the subcritical case (a) and to the supercritical 
case (c) respectively. We represent, for each set of values of {X,Zq), the distance to the 
Yaglom limit, supjgji ... ;^qo} |Pzo(^i = "^1^ < ^o) ~ as a function of the time, and the 
same distance as a function of the logarithm of the survival probability — logP^(,(t < To) 
(i.e. the extinction time scale). By numerical computation, we also obtain that 

(a) A = 0.9: 9{a) = 0.100 and e{a) - x = 0.102. 

(c) A = 1.1: 9{a) = 5.84 x 10"^ and 9{a) - x = 0.103. 

In the case (a), we have 6{a) = 0.100 — x ~ d{c() = 0.102 and we observe that the speed 
of convergence depends on the initial position in a non-trivial way: while the survival 
probability is smaller for the process starting from 10 than for the process starting from 
100, the convergence to the Yaglom limit in the extinction's time scale happens faster 
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in the case Zq = 10. In the case (c), we have 6{a) = 5.84 x 10~^ < X - ^'(a) = 0.103. 
The speed of convergence to the Yaglom hmit in the extinction's time scale depends on 
the initial position: if (□) Zq = 100, then it is almost immediate; if (o) Zq = 10, the 
distance between the conditional distribution and the Yaglom limit is equal to 0.05 when 
the survival probability is around e""'^ ^ 0.61; if (■) Zq = 1, then this distance is equal 
to 0.05 when the survival probability is around e~^'^ ~ 0.091. 

3.2 The Q-process 

Let us now study the marginals of the process conditioned to never be extinct. 

Theorem 9. Assume that we are in the conditions of Theorem^ For any i^, ii, - ■ ■ ,ik & 
E*, any < Si < ■ ■ ■ , < t, the limiting value WiHi^^ 'W^^i^Zg^ — i^, ■ ■ ■ , Zg^ — ^fc|^o ^ 
t) exists. 

Let {Yt,t > 0) be the process starting from ig ^ E* and defined by its finite dimensional 
distributions 

F,,{Yg, =t^,... ,Ys^= tk) = lim Fi,{Zs, =ti,--- ,Zs,= tk\To > t). (14) 
Then Y is a Markov process with values in E* and transition probabilities given by 

It is conservative, and has the unique stationary probability measure {ajnj)j. 

Remark that the stationary probability is absolutely continuous with respect to the QSD, 
but, contrary to intuition, it is not equal to the QSD. 

Proof. Let us denote 6{a) by 6 for simplicity. Let Zq, ii, ■ ■ ■ ,ik & E* and < Si < ■ ■ ■ < 
Sk < t. We introduce the filtration J^s = (^{ZuiU < s). Then 

Pio(^si = ii, - ■ ■ ,Zs^ = ik ; To > t) = (lz,j=n,-,z,^=ife ^io (lro>t|^sj) 

( by Markov property) 
= (Z,, = «i, ■ ■ ■ ,Z,^ = Ik) F,^{To > t - Sk). 



By Theorem |7], 



Thus 



limP,„(Z,, =^i,--- =^fc|To >t) = F,^{Z,^=z,,--- ,Z,^=Zk)'^e''K (15) 



'«0 
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Let us now show that F is a Markov process. We have 



xFi^{Zt-Sk = j) (by Markov property of Z) 



and thus ¥{Yt = j\Y,, = ■ • ■ , F,^ = t^) = ¥i^{Yt_s, = j)- 
By (15) and Theorem [7| we have 



Moreover let us compute the infinitesimal generator L of y from the infinitesimal generator 
L of Z. We have for j ^ i, 



71 ■ 

Lij = limPi,(s) = — Lij. 

s^O TT,- 



For j = i, 



I -Puis) l-e'^'Puis) 
lim ^^ = -lim 



s^O S s 

s^O S 



We thus check that 



Since Lvr = — ^tt, then Xljes* ^i-^'i ~ ^iid thus X]je_B* ~ 0- ^ 

4 QSD for birth and death processes 

We are describing here the dynamics of isolated asexual populations, as for example 
populations of bacteria with cell binary division, in continuous time. Individuals may 
reproduce or die, and there is only one child per birth. The population size dynamics 
will be modeled by a birth and death process in continuous time. The individuals may 
interact, competing (for example) for resources and therefore the individual death's rate 
will depend on the total size of the population. In a first part, we recall and partially 
prove some results on the non-explosion of continuous time birth and death processes. We 
will also recall conditions on the birth and death rates which ensure that the process goes 
to extinction in finite time almost surely. In a second part, we concentrate on the cases 
where the process goes almost surely to zero and we study the existence and uniqueness 
of quasi-stationary distributions. 
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4.1 Birth and death processes 



We consider birth and death processes with rates (Aj)j and that is N-valued pure 

jump Markov processes, whose jumps are +1 or —1, with transitions 



where Aj and yUj, z G N, are non-negative real numbers. 

Knowing that the process is at state i at a certain time, the process will wait for an 
exponential time of parameter Aj before jumping to i + 1 or independently, will wait for 
an exponential time of parameter fii before jumping to i — 1. The total jump rate from 
state i is thus Aj + yUj. We will assume in what follows that Aq = yUo = 0. This condition 
ensures that is an absorbing point, modeling the extinction of the population. Since 
these processes have a main importance in the modeling of biological processes, we study 
in detail their existence and extinction properties, and then their QSDs. 

The most standard examples are the following ones. 

1. The Yule process. For each i G N, Xi = Xi for a positive real number A, and 
/ij = 0. There are no deaths. It's a fission model. 

2. The linear birth and death process, or binary branching process. There exist positive 
numbers A and fi such that Aj = Xi and fii = fii. This model holds if individuals 
reproduce and die independently, with birth rate equal to A and death rate equal 
to /i. 

3. The logistic birth and death process. We assume that every individual in the pop- 
ulation has a constant birth rate A > and a natural death rate fi > 0. Moreover 
the individuals compete to share fixed resources, and each individual j ^ i creates 
a competition pressure on individual i with rate c > 0. Thus, given that the pop- 
ulation's size is i, the individual death rate due to competition is given by c{i — 1) 
and the total death rate is fii = fii + ci{i — 1). 

In the following, we will assume that Aj > and fi^ > for any i G N*. 

We denote by (r„)„ the sequence of the jump times of the process, either births or deaths. 
Let us first see under which conditions on the birth and death rates the process is well 
defined for all time t > 0, i.e. r = lim„ r„ = +oo almost surely. Indeed, if r = lim„ r„ < oo 
with a positive probability, the process would only be defined for t < r on this event. There 
would be an accumulation of jumps near r and the process could increase until infinity 
in finite time. 

Let us give a necessary and sufficient condition ensuring that a birth and death process 
does not explode in finite time. The result is already stated in Anderson [2j, but the 
following proof is actually far much shorter and easier to follow. 

Theorem 10. The birth and death process does not explode in finite time, almost surely, 
if and only if rn = +oo, where 



2 — )■ 2 + 1 with rate Xi 
2 — )■ i — 1 with rate fii 




n-l 
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Proof. 1) Let us more generally consider a pure jump Markov process {Xt,t > 0) with 
values in N, and generator {Lij,i,j G N). We set = —La. Let (r„)„ be the sequence of 
jump times of the process and {Un)n the sequence of inter-times defined by 

Un = rn- Tn-l, V?2 > 1; Tq = 0, Uq = 0. 

We also set Too = lim„_>oo Tn G [0, + oo]. The process does not explode in finite time 
almost surely (and is well defined for all time t G if and only if for each i G N 

Pi (Too < oo) = 0. 

Let us show that this property is equivalent to the fact that the unique non-negative and 
bounded solution x = (a;i)ieN of Lx = x is the null solution. 



For any i, we set hf*^ = 1 and, for n G N*, hi"' = Ej(exp(— X]fc=i Uk))- For any n G N, 
we have 



An) 



h^^'^ = J2^hf E,(exp(-f/0). 



Indeed, the property is true for n = since Xlj^j ~^ ~ ^- Moreover, by conditioning with 
respect to Ui and using the strong Markov property, we get 



n+l 



exp(- Uk 



k=l 



exp(-f/i) Exy^ exp(-^[/fc 



(16) 



k=l 



since the jump times of the f/i-translated process are the r„ — Ui,n G N*. We have 



J2HXu,=j) EJexp(-^t/fc 



k=l 



k=l 



5^-^E,(exp(-^^,)), 



k=l 



since Pj(Xc7^ = j) = — . By independence of Ui and Xj/^, we deduce from (16) that 



n+l 



E, exp(- J^Uk)] =J2^^j\ exp(- ^ U^) E^ (exp(-f/i)) . 



k=l 



k=l 



As 



it turns out that 



E,(exp(-f/i)) = / g.e-^'^e-^s ^ , 



(n+l) 



+ qi 



Let (xj)j be a nonnegative solution of Lx = x bounded by 1, then Xj 

LiiXi + -^ij-^j ~ Qi^i "I" -^ij-^ji that 



(17) 
(18) 
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Since h\ = 1 > Xi > and > for all i,j G E, we deduce by iteration from (17) 
and (pl that h^^ >Xi>0, for any neN. 



Let us in the other hand define for any j the quantity Zj = E,j{e ^°°). Using t^q = lini„ r„, 
and r„ = Y12=i we deduce by monotone convergence that zj = lim„ h^^\ 

If the process does not explode a.s., then t^o = oo a.s., and lim„ /i^-"''' = Zi = 0. Since 
^i"'* > a^j > 0, we deduce that Xj = 0. It turns out that the unique nonnegative bounded 
solution of Lx = x is zero. 

If the process explodes with positive probability, then there exists i such that Pj(roo < 

L, 



oo) > 0. Making n tend to infinity in (17), we get Zi = Ylj^i 1+^ ^j- Since > 0, z is a 



positive and bounded solution of Lz = z. 

2) Let us now apply this result to the birth and death process with Aq = /io = 0. Then 
for i > 1, = Aj, Lii_i = fii, Li^i = — (Aj + /ij). The equation Lx = x is given by 

Xq = and for all n > 1 by 

An-^n+l (An ~l~ fJ'n)'^n ~l~ ^'U'^n—l ■^n- 

Thus, if we set A„ = x„ — x„_i, we have Ai = Xi and for n > 1, A„+i = A^ + -^Xn- 
Let us remark that, for any n, A„ > and the sequence {Xn)n is nondecreasing. If Xi = 0, 
the solution is zero. If not, we get by induction 



A„+i = —Xn + 2_^Y~ ~\ ^ + ^; — a^i- 

Letting 



^ n— 1 

^ _ -1- _l_ /^fc+l ■ • ■ /^n _|_ A'l • ■ ■ A'n. 

An j^^-^ AfcAfc+i ■ ■ ■ An Ai ■ ■ ■ An 

we deduce that rn Xi < An+i < Tn Xn- Then 



Xi(l + ri H Tn) < Xn+l < Xi JJ(1 + Tfc). 

The boundedness of the sequence (x„)„ is thus equivalent to the convergence of the series 
Ekrk. □ 

Corollary 11. Let us consider a BD-process with birth rates {Xi)i. If there exists a 
constant A > such that 

K < Ai, \fi > 1, 
then the process is well defined on IR+. 

The proof is immediate. It turns out that the linear BD-processes and the logistic pro- 
cesses are well defined on IR+. 

Let us now recall under which assumption a BD-process goes to extinction almost surely. 
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Proposition 12. The BD-process goes almost- surely to extinction if and only if 



oo 



Proof. Let us introduce 

:= P(£'xtmct2on|Zo = i) = Pi(To < oo), 

which is the probability to attain in finite time, starting from i. As before Tq denotes 
the extinction time and T/ the hitting time of any /. The Markov property yields the 
induction formula 

\i Ui+i - (Ai + /ii) Ui + jji = 0, Vi > 1. 

To resolve this equation, we firstly assume that the rates Xi, fit are nonzero until some 
fixed level / such that Xj = fij = 0. We set for each i, u'p := Pi(To < Tj). Thus 
Ui = lim/^oo ^^1^^ • Defining Uj := Ylk=i ^x~^^ easy computation shows that for i G 
{I,-- - ,/-!}, 

7-1 



k=i 



Ai ■ ■ ■ A 



k 



In particular, -u^^^ = y+jT^- Hence, either (f//)/ tends to infinity when / — )■ oo and any 
extinction probability Ui is equal to 1 or (?77)7 converges to a finite limit Uoo and for i > 1, 



oo 



(i+f/oo)-EF^<i- 

fe=J 



□ 



Corollary 13. 1. The linear BD-process with rates Xi and fii goes almost surely to 
extinction if and only if X < fi. 

2. The logistic BD-process goes almost surely to extinction. 

Proof. 1) If A < fi, i.e. when the process is sub-critical or critical, we obtain Uj > I — 1 
for any / > 1. Then {Ui)i goes to infinity when / — )■ oo and the process goes to extinction 
with probability 1. Conversely, ii X > fi, the sequence {Ui)i converges to and an 
easy computation gives Ui = {X/fiy. 

2) Here we have 

Xi = Xi ; fii = fii + ci{i — 1). (20) 



It is easy to check that (19) is satisfied. □ 
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4.2 Quasi-stationary distributions for birth and death processes 

We consider a BD-process (Zt) with almost sure extinction. A probability measure a on 
N* is given by a sequence {<yj)j>i of non-negative numbers such that Ylj>i O-j = ^■ 
Our first result is a necessary and sufficient condition for such a sequence {o.j)j>i to be a 
QSD for Z. Thereafter we will study the set of sequences which fulfill this condition (we 
refer the reader to van Doorn [59] for more details). 

Theorem 14. The sequence {aj)j>i is a QSD if and only if 

1. ttj > 0, Vj > 1 and J2j>i^j — 1- 

2. yj> 1, 

Aj_iaj_i — {Xj + fJ.j)aj + /ij+iOj+i = — /iiaia^ ; 

— (Ai + + /i2tt2 = — /iiOi- (21) 



The next result follows immediately. 

Corollary 15. Let us define inductively the sequence of polynomials {Hn{x))n as follows: 
Hi{x) = 1 for a// X G M and for n > 2, 

Xn Hn+l{x) = (A„ + fln-x) Hn{x) - f^n-l Hn-l{x) ] 

Ai H2{x) = Xi+iii-x. (22) 
Then, any quasi- stationary distribution {aj)j satisfies for all j > 1, 

Oj = ai TTj Hj{jj,iai), 

where 

Ai ■ ■ ■ A„_i 

TTi = 1 ; nn = . (23) 

■ ■ ■ /in 



Proof of Theorem 14 By Proposition |4] and for a QSD a, there exists 6 > such that 

aL = —9 a, 

where L is the infinitesimal generator of Z restricted to N*. Taking the j^^ component of 
this equation, we get 

Aj_iaj_i — {Xj + jJ.j)cij + /ij+iaj+i = —6aj, Vj > 2 
— (Ai + fii)ai + /i2«2 = —9ai- 

Summing over j > 1, we get after re-indexing 



We deduce that 9 = /iiai, which concludes the proof of Theorem 14 



□ 
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The study of the polynomials [Hn) has been detailed in Van Doom [59]. In particular it 
is shown that there exists a non-negative number such that 

X < ^1 Hnix) > 0, Vn > 1. 



By Corollary 15, aj = ai Hj{niai). Since for any j, aj > 0, we have Hj{jj,iai) > for 
all j > 1 and then 

< /iifti < ^1. 

We can immediately deduce from this property that if = 0, then there is no quasi- 
stationary distribution. 

To go further, one has to study more carefully the spectral properties of the semi-group 
(Pt) and the polynomials (if„)„, as it has been done in ^39], [30j and [59j. From these 
papers, the polynomials {Hn)n are shown to be orthogonal with respect to the spectral 
measure of {Pt}- In addition, it yields a tractable necessary and sufficient condition for 
the existence of QSD based on the birth and death rates. The series {S) with general 
term 

^ oo 

Sn = \ / VTj 

i=n+l 

plays a crucial role. Remark that {S) converges if and only if 



< +00. 



Theorem 16. (159], Theorems 3.2 and 4.1). We have the convergence 

lim ¥i{Zt = j\To >t) = — Tij ^1 H,{^i). 

In particular, we obtain 

^1= lim/iiPi(Zt = l|To>t) (24) 

1. If C,i = 0, there is no QSD. 

2. If (S) converges, then > and the Yaglom limit is the unique QSD. 

3. If (S) diverges and ^ 0, then there is a continuum of QSD, given by the one 
parameter family {&j{x))o<x<$,i 

aj{x) = — TTj X Hj{x). 



Remark 4. 1. Formula (24) and the approximation method described in Section ro^ 



allow us to deduce a simulation algorithm to get ^i. 
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2. Cases with more general birth and death processes have also been studied recently. 
Let us mention the infinite dimensional state space setting of Collet, Martinez, 
Meleard and San Martin [15], where each individual has a type in a continuous 
state space which infiuences its birth and death rates, and mutation on the type can 
occur. The authors give sufficient and quite general conditions for the existence and 
uniqueness of a QSD. We also refer the reader to the recent work of van Doorn j60j . 
where a transition to the state may occur from any state. The author provides 
sufficient conditions for the existence of QSDs. 



Let us now develop some examples. 



The linear case. We assume Aj = Az ; = fxi and A < /i. In that case, the BD-process 
is a branching process, where each individual reproduces with rate A and dies with rate 
/i. A straightforward computation shows that the series (S) diverges. 

Setting fs : k ^ s^, we get by the Kolmogorov forward equation, 

^^-^ = ^^PtfsiO) - (A + /i)P,/.(l) + AP,/,(2). 

But the branching property of the process implies Ptfs{2) = (Pt/s(l))^, while /s(0) = 1 
so that 



dt 

2(1 - s)(2m - 1) 



Setting m = 2yt-, we deduce that for s < 1 



Ptfs{l) = l' 



{ms + m- 2)e-(^+^)(™-i)* + (1 - s) 



In particular, we deduce that the generating function Ft : s ^ K{s^^\Zt > 0) of Zt 
conditioned to > converges when t goes to infinity: 

P,/,(l) - P,/o(l) {\-fi)s 

*^ ^ l-P*/o(l) As-/i- 

We deduce that the Yaglom limit of Z does not exist if A = and is given by the geometric 
distribution with parameter ^ if A < /i: 

xV-'r^ A 



An easy computation yields = /i — A, since by (24), ai = ^. But the series (S) diverges 

so that for A < yU, > and there is an infinite number of QSD. If A = /x, i^i = and 
there is no QSD. 



The logistic case. We assume Aj = Ai ; /ij = /li + d(i — 1). Because of the quadratic 
term, the branching property is lost and we can not compute the Yaglom limit as above. 
Therefore, we have no other choice than to study the convergence of the series (S). 
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We have 

2^ - ^Ici i\~^\c) (n + p+l)\ 

i=n+l i=n+l ^ ' p=0 ^ ^ \ < 1^ < J 

- [cj (n + l)!^U; " \c) {n + l)\ 
since ^^gg^ < ^. Thus as ^ < C (f n! , we get 

C 1 A 



5Z 



< 



A„7r„ c + 1) 

i=n+l 



Hence the series converges. Thus the Yaglom hmit exists and is the unique quasi- 
stationary distribution. 

One can obtain substantial quahtative information by looking closer to the jump rates 
of the process. For instance, (A — /i)/c is a key value for the process. Indeed, given a 
population size i, the expectation of the next step is equal to j[tM^^^+ci(^-\\^ • Then the sign 
of this expectation depends on the position of z — 1 with respect to (A — /i)/c: 

• If i < (A — /i)/c + 1, then the expectation of the next step will be positive. 

• If 2 = (A — /i)/c + 1, then it will be 0. 

• If 2 > (A — /i)/c + 1, then it will be negative. 

We deduce that the region around (A — yu)/c is stable: it plays the role of a typical size for 
the population and we expect that the mass of the Yaglom limit is concentrated around 
it. The value (A — /i)/c is called (by the biologists) the charge capacity of the logistic BD- 
process with parameters A, n and c. In the next section, we will consider large population 
processes, which means logistic BD-processes with large charge capacity. 

Example 3. We develop now a numerical illustration of the logistic BD-process case. 
Across the whole example, the value of the charge capacity is fixed, arbitrarily chosen 
equal to 9. 

In order to illustrate the concept of charge capacity, we represent in Figure |4] a random 
path of a logistic birth and death process with initial size Zq = 1 and with parameters 
A = 10, /i = 1 and c = 1. We observe that the process remains for long times in a region 
around the charge capacity. Moreover, we remark that the process remains mainly below 
the charge capacity; this is because the jumps rate are higher in the upper region, so that 
it is less stable than the region below the charge capacity. 

Let us now compare the Yaglom limits (numerically computed using the approximation 
method presented in Section [gJ of two different logistic BD processes whose charge capac- 
ities are equal to 9 (see Figure [s]): 



(a) Z^"'\ whose parameters are A = 10, /i = 1 and c = 1, 

(b) Z^''\ whose parameters are A = 10, /x = 7 and c = 1/3. 
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Figure 5: Example 3. The Yagloni limits of two logistic birth and death processes with the same charge 
capacity = 9: (a) A = 10, = 1 and c = I; (h) X — 10, n = 7 and c — 1/3. 
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Z0=100 
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Z0=10 









Figure 6: Example 3. Evolution of the distance between the conditioned distribution and the Yaglom 
limits of two logistic birth and death processes with the same charge capacity — 9: (a) A = 10, /i = 1 
and c = 1; (b) A = 10, /i = 7 and c = 1/3. 



We observe that the Yaglom hmits of Z^") and Z^''^ are supported by a region which is 
around the charge capacity. We also remark that the Yaglom limit of the process 
has a more flat shape than the Yaglom limit of This is because the competition 

parameter of Z^^^ is small in comparison with the birth and death parameters, so that the 
drift toward the charge capacity is small too, both above and below the charge capacity. 

We compute now the distance between the conditioned distribution and the Yaglom limit 
for the two processes Z*-") and Z'-^-' for different values of the initial state, namely Zq = 1, 
Zq = 10 and Zq = 100. The numerical results are represented in Figure |6} We observe 
a strong dependence between the speed of convergence and the initial position of the 
processes. In the case of Z^"'\ it only takes a very short time to the process starting from 
100 to reach the charge capacity, because the competition parameter is relatively high 
and so is the drift downward the charge capacity. On the contrary, in the case of Z^^\ it 
takes a longer time for the process to come back from 100 to the charge capacity, so that 
the speed of convergence to the Yaglom limit is slow. In both cases, the convergence to 
the Yaglom limit happens very fast when starting from the value 10, because it is near 
the charge capacity. 
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5 The logistic Feller diffusion process 



5.1 A large population model 

We are now rescaling logistic birth and death processes with a parameter K modehng 
a large population with small individuals, i.e. a population with a large initial size of 
order K and a large charge capacity assumption. The individual's weights (or biomasses) 
are assumed to be equal to and we study the limiting behavior of the total biomass 

process {^,t > 0) when K tends to infinity, being the population's size at time t. 
In what follows. A, fi and c are fixed positive constants. 



In Subsection |5.1.1[ individual birth and death rates are assumed to be constant and 
the competition rate depends linearly on the individual biomass In Subsection 



5.1.2 



we investigate the qualitative differences of evolutionary dynamics across populations 
with allometric demographics: life-lengths and reproduction times are assumed to be 
proportional to the individual's weights. 

In both cases, the charge capacity of (Z^) will be (A — fi)K/c. 
5.1.1 Convergence to the logistic equation 

Given a parameter K scaling the population's size, we consider the logistic BD-process 
Z^ with birth, death and competition parameters A, /i and c/K respectively. We assume 
that the initial value of Z^ is of order K, in the sense that there exists a non-negative 
real random variable Xq such that 

Z^ 

> Xq with E(Xo^) < +00. 

We consider the total biomass process defined by X^ = Z^^ /K for all i^' > 1 and are 
interested in the limit of X^ when K — )■ oo. The transitions of the process {X^,t > 0) 
are the following ones: 

— — with rate Xi = XK — ; 

K K K' 

^ ^ withrate /i2 + |,2(z-l) = ir-^(^/i + c(-^--^: 

Theorem 17. Assume that Xq is a positive number xq. Then, the process {X^,t > 0) 
converges in law m D([0,T], IR+) to the unique continuous (in time) deterministic function 
solution of 

x{t) = xo -l- / (A — /i — cx{s))x{s)ds. 
Jo 

Remark 5. The function x is thus solution of the ordinary differential equation 

X = (X — n)x — cx^ ; x(0) = Xq, (25) 
called logistic equation. This equation has been historically introduced as the first macro- 



scopic model describing populations regulated by competition (cf. iM])- In Theorem 17 
above, it is obtained as the limit of properly scaled stochastic jump models. 
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The function x solution of (25) hits in finite time if A < /i, while it remains positive 
forever if A > yU, converging in the long term to its charge capacity Thus at this 

scale extinction does not happen. 



Proof of Theorem 11. The Markov process {X^ , t > 0) is well defined and its infinitesimal 
generator is given, for any measurable and bounded function 0, by 

Lk0(x) = \Kx + ^)- + K{^ix + cx{x - ^)) (^^{x - ^) - ■ (26) 

Hence, by Dynkin's theorem ([21] Prop. IV-1.7), the process 

</.(Xf ) - 0(Xo^) - fLK<l>{X^)ds 
Jo 



(27) 



defines a local martingale, and a martingale, as soon as each term in (27) is integrable. 
In particular, taking (f){x) = x gives that (Xf ,t > 0) is a semimartingale and there exists 
a local martingale such that 



ds. 



(28) 



Since Xq is deterministic and using a localization argument, we deduce that E(supj<'p(X/^^) ) < 



oo. Moreover, taking (f){x) = x^ applied to (27), and comparing with Ito's formula ap- 
plied to {X^Y prove that {M^) is a square-integrable martingale with quadratic variation 
process 

{M% = i /7a + /X + c (xf - ^\ \ Xfds. (29) 



K 



K 



Let us now study the convergence in law of the sequence {X^), when K tends to in- 
finity. For any the law of X^ is a probability measure on the trajectory space 
Dj, = D([0,T], ]R+), namely the Skorohod space of left-limited and right-continuous func- 
tions from [0,T] into IR+, endowed with the Skorohod topology. This topology makes 
a Polish state, that is a metrizable complete and separable space, which is not true if 
is endowed with the uniform topology. See Billingsley |9j for details. 



ness of the solution of (25) is immediate. 



The proof of Theorem 17 is obtained by a compactness-uniqueness argument. The unique- 



By a natural coupling, one may bound the birth and death process X^ stochastically 
from above by the Yule process started from Xq, which jumps from x to x + at the 
same birth times than X^ . One easily shows that sup^ E(supi<j'(Fj^)^) < oo and thus 

supE(sup(Xf )^) < oo. 



K 



t<T 



From this uniform estimate, we deduce the uniform tightness of the laws of X (as 
probability measures on D^), using the Aldous criterion (cf. Aldous [I], Joffe-Metivier 
[37J). Then, by Prokhorov's Theorem, the compactness of the laws of {X^) follows. Since 
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is continuous on D^, 



< ^ and the function x h- )■ sup^^^ \xt — Xt- 
each limiting value (in law) of the sequence {X^) will be a pathwise continuous process. 
In addition using (29) and ( |5.1.1 ), it can be shown that \mvK^^^{{M^)t) = 0. Then, 
the random fluctuations disappear when K tends to infinity and the limiting values are 
deterministic functions. Now it remains to show that these limiting values are solutions 
of (25), which can be done similarly to the proof of Theorem 18 (4) stated below. □ 



5.1.2 The logistic Feller diffusion process 

In this section, we study the logistic BD-processes with birth and death rates given 
by 7 + A and 'y K + fi respectively. Here 7, A and /i are still positive constants. We 
assume that the competition parameter is given by c/K, so that the charge capacity of 
is still (A - fi)K/c. 

Remark 6. This BD-process {Z^)t can also he interpreted as a time-rescaled BD-process 
whose birth, death and competition parameters are given by 'j + \/K , 7 + yu/i^ and 
c/K'^ respectively, that is a critical BD-process with small pertubations. 



We are considering as in Section 5.1.1 the sequence of processes defined for alH > 
by = f . 

The transitions of the process {X^) are given by 

i i + l 

K 

i 

K 



K 

I - 1 



K 



with rate '-yKi + \i 



with rate '-)Ki + /iz H i{i 

K 



(30) 



Formula (28) giving the semi-martingale decomposition of X^ will stay true in this case 
with another square integrable martingale part such that 



1 l\2^K + A + /i + c (^Xf - 1^ X^ds. 



One immediately observes that the expectation of this quantity does not tend to zero as 
K tends to infinity. Hence the fluctuations will not disappear at infinity and the limit 
will stay random. Let us now state the convergence theorem. 

Theorem 18. Assume 7, c. A, /i > and X > fi. 

i) Assume that the sequence {Xq )k converges in law to Xq with IE(Xq) < oc. Then 
the sequence of processes {X^)k with transitions (30) converges in law in V{Vit) to the 
continuous process X , defined as the unique solution of the stochastic differential equation 

dXt = ,/2^tdBt + ((A - fi)Xt - cXf) dt- Xo G]0, + oo[, (31) 

where {Bt)te[o,+cx,[ is a standard Brownian motion, 
a) Let us introduce for each y > the stopping time 

Ty = M{teR+,Xt = y}. (32) 

For any a; > 0, we get 

p,(ro < 00) = 1. 
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When c = 0, Equation (31) defines the super-critical Feller diffusion process and will 
explode with positive probability. In the general case where c 7^ 0, it defines the so-called 
logistic Feller diffusion process following the terminology introduced by Etheridge |2QJ and 
Lambert [42j . Let us remark that the quadratic term driven by c regulates the population 
size, which fluctuates until it attains the absorbing point 0. Theorem 18 shows that 
an accumulation of a large amount of birth and death events may create stochasticity, 
often called by biologists ecological drift or demographic stochasticity. Contrarily to the 
previous case (Theorem 17), the limiting process suffers extinction almost surely. 



Proof. As for Theorem [17], the proof is based on a uniqueness-compactness argument. 

(1) The uniqueness of the solution of (31) follows from a general existence and pathwise 
uniqueness result in Ikeda-Watanabe ^36j Section IV-3 or Karatzas-Shreve [3§J. For a 
stochastic differential equation 

dXt = cT{Xt)dBt + b{Xt)dt, 

with a and b smooth enough, the existence and pathwise uniqueness are determined 
thanks to the following scale functions: for a; > 0, 



Q{x) 



oQ{y) 



dy ; A(a;) 



-'^^'Uz ] dy. 



(33) 



More precisely, it is proved that (i) Vx > 0, Px(7o < Too) = 1 and (ii) A(+oo) = 
+00 ; /t(0+) < +00 are equivalent. In that case, pathwise uniqueness follows, and then 

uniqueness in law. 

In our situation, the coefficients are given by 

cr(x) = \/2rfx ; b{x) = (A — fi)x — cx"^, 

so that the functions A and n satisfy (ii). Thus the SDE (31) has a pathwise unique 
solution which reaches in finite time almost surely. 

(2) Let us assume that IE(A'q) < 00 and prove that supj^ E(sup^<2-(X/'')^) < 00. The 
infinitesimal generator of X^ is given by 

1 



LK(j){x) = {-fKx + Xx)K { (j){x + —) - (pix) 



+ {'yKx + fxx + cx{x — -^))K I 0(x 



K' 



1. 
K' 



(t){x] 



(34) 



With 0(x) 
+ 



X 



we obtain that 



^0^ + + l iK'Xf [xf + 1) - (xf - 1) - [X^f 



ds 



XKX^ 



Xf + 



K 



+ / {ii + c{Xf -l))KX^ 



X 



K 



ds 
1 

K 



K\3 



ds. 



35 



where is a local martingale. Using a standard localization argument and 



(xf + ^)'-(xf-l)'-(Xfr = 6 



X. 



K 



we get 



Jo 

where C is independent of K. Gronwall's lemma yields 

supsupE(|Xf H < oo. 



t<T K 



(35) 



Now, thanks to (|35j) and Doob's inequality, we deduce from the semi-martingale decom- 
position of {X^Y (obtained using (j){x) = x^), that 



supE(sup |Xj 

K t<T 



K\2\ 



< OO. 



(36) 



(3) As previously, the uniform tightness of the laws of {X^) is obtained from (36) and the 



Aldous criterion [Ij. Therefore, the sequence of laws is relatively compact and it remains 
to characterize its limit values. 



(4) As in the proof of Theorem 17, we remark that the limiting values only charge the 
set of continuous trajectories, since sup^^j. | AX/^| < j^. Let Q G P(C([0,T], M+)) be a 
limiting value of the sequence of laws of the processes X^^. We will identify Q as the 



(unique) law of the solution of (31) and the convergence will be proved. Let us denote 
Ct = C([0,r],M+) and define, for (j) e and t > 0, the function 

X ^ 0(XJ - 0(Xo) - /" (7X,0"(X,) + ((A-/x)X,-cX2)0'(X,))rfs, 

Jo 

which is continuous Q-a.s.. Let us show first that the process {tpt{X))t is a Q-martingale. 
For X G M+, we define 

L(j){x) = 7X0" (x) + ((A — fi)x — cx^)0'(x). 



Using Taylor's expansion, we immediately get (with defined in (34)) 

'jK'^ X 
+X K X 



(x)-L0(x)| 



0(x + -i) + 0(x - -^) - 20(x) - ^J\x) 



(x + -^) - (p{x) - ^<P'{x) 



< 



+K (/ix + cx(x — 1)) 
C 



1 1 

(x - — ) - 0(x) + j^(p'{x) 



K 



(37) 
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where C doesn't depend on x and K. By ([36]), we deduce that E (^|Li^0(Xf ) - 
tends to as if tends to infinity, uniformly for t G [0,T]. 

For Si < ■ ■ ■ < Sfc < s < t, for Qi, - ■ ■ ,gk G C*6, we introduce the function H defined on 
the path space by 

HiX) = g,{X,J ■ --gkiXj - ^,(X)) . 

Let us show now that 

Eq{H{X)) = 0, (38) 
which will imply that {tlJt{X))t is a Q-martingale. 

By construction, ^f(X^) = 0(X/^) — 0(Xo) — /J LK(p{Xf)ds defines a martingale, then 

E ■ ■ ■ g.iX^l) (^f (X^) - V^f (X^))] = 0. 

In another way, this quantity is equal to 

E [giiXfj . . . g,{Xfj (^f (X^) - V^f (X^) - MX'') + ^.(X^))] 

+E [g^iXfj . . . g,{Xfj (^,(X^) - ^,(X^)) - g,{X,,) ■ ■ ■ g,{X,,) (MX) - ^.(X))] 

+E [^i(X,J ■ ■ ■ ^,(X,J (V^i(X) - ^,(X))] . 



and tends 



The first term is equal to E [^i(Xj^) ■ ■ ■ gkiXfJ £ (^L;^0(Xf ) - L0(Xf ) j (is 
to by (g and (|37|. 

The second term is equal to K{H{X^) — H{X)). The function X i-> H{X) is continuous 
and since H{X) < C (l + f^{l + X^)du}j, it is also uniformly integrable by (Q. This 
leads the second term to tend to as if tends to infinity. 

Therefore, it turns out that (38) is fulfilled and the process iptiX) = (f){Xt) — 0(Xo) — 
L(j){Xs)ds is a Q-martingale. 

By (Q and taking (f){x) = x leads to X^ = Xq + + /J((A - fi)Xs - cX'^)ds, where M 
is a martingale. Taking (f){x) = x"^ on the one hand and applying Ito's formula for X^ on 
the other hand allow us to identify 

{M)t= f^'jX.ds. 



By the representation theorem proved in |38| Theorem III-4.2 or in [36], there exists a 
Brownian motion B such that 



Jo 



That concludes the proof. 



□ 
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5.2 QSD for logistic Feller diffusion processes 
5.2.1 Statement of the results 

We are now interested in studying the quasi-stationarity for the logistic Feller diffusion 
process solution of the equation 



dZt = VZtdBt + (rZt - cZt)dt, Z^ > 0, 

where the Brownian motion B and the initial state Zq are given, and r and c are assumed 
to be positive. (We have assumed that 7 = 1/2). The results and proofs that are presented 



in Section 5.2 have been obtained by Cattiaux, Collet, Lambert, Martinez, Meleard and 
San Martin 



Let us firstly state the main theorem of this part. 

Theorem 19. Assume that Zq, r and c are positive. Then the Yaglom limit of the process 
Z exists and is a QLD for Z starting from any initial distribution. As a consequence, it 
is the unique QSD of Z . 

Remark 7. 1) The theory studying the quasi-stationary distributions for one-dimensional 
diffusion processes started with Mandl |45J and has been developed by many authors. See 
in particular [16|, |17], [37], [40]. Nevertheless in most of the papers, the diffusion and drift 



coefficients are regular and the "Mandl's condition" k{+oo) = 00 (see (33)) is assumed. 
This condition is not satisfied in our case because of the degeneracy of the diffusion and 
the unboundedness of the drift coefficient. 



2) Theorem 19 differs from the results obtained in case of drifts going slower to infinity. 
For example, Lambert proves that if c = and r < 0, then either r = and there 
is no QSD, or r < and there is an infinite number of QSD. Lladser and San Martin 
[44J show that in the case of the Ornstein-Uhlenbeck process dYt = dBt — Ytdt, killed at 
0, there is also a continuum of QSD. In the logistic Feller diffusion situation as in the 
logistic BD-process, the uniqueness comes from the quadratic term cX^ induced by the 
ecological constraints. 

3) We have seen that the rescaled charge capacity of the logistic birth and death process 
converges to the charge capacity of the logistic Feller diffusion. However, whether the 
rescaled Yaglom limit of the logistic birth and death process converges to the Yaglom 
limit of the logistic Feller diffusion process remains an open problem. 



In order to prove Theorem [T9| we firstly make a change of variable and introduce the 
process {Xt,t > 0) defined by Xf = 2a/Z^. Of course, X is still absorbed at and QSDs 
for Z will be easily deduced from QSDs for X. From now on, we focus on the process 
(Xt). 

An elementary computation using Ito's formula shows that (Xt) is the Kolmogorov diffu- 
sion process defined by 

dXt = dBt - qiXt)dt, (39) 

with 

1 rx cx^ 

^(^) = 2^-y + x- 
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Mention that the function q is continuous on but explodes at as and at infinity as 
The strong (cubic) downward drift at infinity will force the process to live essentially 
in compact sets. That will provide the uniqueness of the QSD, as seen below. 

We introduce the measure fi, defined by 

where Q is given by 

Q{y) = J' 2q{z)dz = \ny + "^{1 - y') + ^{y' - 1). (40) 

In particular —Q/2 is a potential of the drift —q. The following result clearly implies 
Theorem [l9l 



Theorem 20. 1 13^ Assume that Xo, r and c are positive. Then the Yaglom limit a of 
the process X exists. 

Moreover, there exists a positive function rji G lJ{dfi) such that 

1. 

c^idx) = r ^ ,f gfa) , dx, (41) 

2. Vx e M;, lim^^o, e^(")*Px.(To > t) = r]i{x), 

3. there exists x > such that, Vx G W^, 

lim e-^'^-^^"))* |Px. {Xt G A\To > t) - a{A)\ < +00. 

4. the QSD a attracts all initial distribution, which means that a is a QLD for X 
starting from any initial distribution. 



The proof of Theorem 20 will be decomposed in the next subsections. 



5.2.2 Spectral theory for the killed semi-group 



As previously we are interested in the semi-group of the killed process, that is, for any 
X > 0, t > and any / G Cb(R\), 

PJ{x)=E,{f{Xt)h^To), (42) 
with the associated infinitesimal generator given for G C^((0, + 00)) by 

= ^0" - 

We are led to develop a spectral theory for this generator in L^(yu). Though the unity 
function 1 does not belong to L^(/i), this space is the good functional space in which to 
work. The key point we firstly show is that, starting from x > 0, the law of the killed 
process at time t is absolutely continuous with respect to fi with a density belonging to 
L^(/i). The first step of the proof is a Girsanov Theorem. 
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Proposition 21. For any bounded Borel function F defined on Q = C([0,t],]R!j_) it holds 



E, [F{u)lt<Toico)] =E^ 



/I 1 1 /■ 

F{uj)lt<ToHexp i-Q{x) - -Q{uJt) - - (g^ - q' 



){ujs)ds 



where K^"^ denotes the expectation with respect to the Wiener measure starting from x 
and u the current point in Q. 

Proof. It is enough to show the result for non-negative and bounded functions F. Let 
e G (0,1) and = A Ti/^. Let us choose some ip^ which is a non-negative C°° function 
with compact support included in ]e/2,2/e[ such that '4^e{u) = 1 if e < m < 1/e. For all 
X such that e < x < 1/ e the law of the diffusion (39) coincides up to with the law of 
a similar diffusion process obtained by replacing q with the cutoff function q^ = qip^. 
For the latter we may apply Novikov criterion (cf. [53j p. 332), ensuring that the law of 
is given via Girsanov's formula. Hence 



E, [F{u)lt<r,^^)] =E^ 

= E^ 
= E^ 



F{u) lt<r,{u^) exp (^J^ -q 



-qe{ujs)dus 



F{uj) lt<r,{u,) exp 



-q{us)dujs - 



{qef{us)ds 



q^{us)ds 



(1 1 1 /■ 

^(^)li<r.H exp \-Q{x) - -Q{ujt) - - (g^ - q' 



){ujs)ds 



integrating by parts the stochastic integral. But li<T-^ is non-decreasing in e and converges 
almost surely to lt<To both for W^. and for P^, (since Fx{Tq < oo) = 1)). Indeed, almost 
surely. 



limX^^ 



limX^^ 



lim e = 

£->0 



SO that lim£_j.o > Tq. But < Tq yielding the equality. It remains to use Lebesgue 
monotone convergence theorem to finish the proof. □ 

Theorem 22. For all x > and all t > there exists a density function r{t,x,.) that 
satisfies 

r+oo 

E4f{Xt) lt<ro] = / f{y) r{t,x,y) fi{dy) 
Jo 

for all bounded Borel function f . In addition, for all t > and all x > 0, 

X) 

r\t,x,y) i^idy) < (l/27rt)^ e^* e«(") , 



where 

Proof. Define 



C = -m{{q\y)-q'{y)) < +oo. 

y>0 



fl 1 1 /"* 

lt<ToM exp ( - Q{uJo) - 2 Q{^t) ~ 2 Jo ^ ^' 



- q'){ujs)ds 
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Denote by e'^'^*'^'?') = (27rt)- "2 exp I — J!' ) the density at time t of the Brownian motion 
starting from x. According to Proposition [2T| we have 

E.(/(X,)WJ = E*^ (/(a;,) E*^(G|a;,)) 

r+oo 

Jo 

r+oo 

= fiv) E*^ (Glut = y) e-''(*'^'^)+«(^) /i(rfy) , 

^0 

because Fy^'^lGlut = y) = if ?/ < 0. In other words, the law of Xt restricted to non 
extinction has a density with respect to fi given by 

r{t,x,y) = EW-(G|wt = y) e-^^'^'y^+'^^y^ . 

Hence 

+00 /• + 00 2 



r\t,x,y) fi{dy) = / (E*-(G|a;i = e-^(*'^'^')+«(^)) 

= E^- (EW-(G'|u;i))') 
< E^" (e"''^*'''''^*^^'^'^'^'^ E^"(G^|u;t)) 

where we have used Cauchy-Schwarz's inequahty. Since e~'"^^'^''^ < (l/2nt)^, the proof is 
completed. □ 



Thanks to Theorem 22 , we can show, using the theory of Dirichlet forms (cf. Fukushima's 



book |25| ) that the infinitesimal generator L of X, defined by (5.2.2), can be extended 
to the generator of a continuous symmetric semi-group of contractions of L^(yu) de- 
noted by {Pt)t>o- In all what follows, and for f\g E L^(yu), we will denote {f,g)^ = 
J^^ f{x)g{x)n{dx). The symmetry of Pt means that {Ptf,g)^, = {f,Ptg)f,. 
In Cattiaux et al. |13] . the following spectral theorem in L^(yu) is proved. 

Theorem 23. The operator —L has a purely discrete spectrum < Ai < A2 < ■ ■ ■ • 

Furthermore each Aj (i G W) is associated with a unique (up to a multiplicative constant) 
eigenfunction rji of class C^((0,oo)), which satisfies the ODE 

^Vi - QVi = -^iVi- (43) 

The sequence {r]i)i>i is an orthonormal basis o/L^(/i) and rji^x) > for all x > 0. 
In addition, for each i, rji G 'L^{^). 
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The proof of this theorem is based on a relation between the Fokker-Planck operator L 
and a Schrodinger operator. Indeed, let us set for g G h'^{dx), 

P,g = e'^'^ P,{g e^l\ 

Pt is a strongly semi-group on L^((ix) with generator defined for g G C^((0, + oo)) by 

The spectral theory for such Schrodinger operator with potential ~^ on the line (or 
the half-line) is well known (see for example the book of Berezin-Shubin [7j), but the 



potential 



does not belong to as generally assumed. Nevertheless, in our case 



inf(g^ — q') > — oo, which ensures the compactness of the operators L and Pt. 



The following corollary of Theorem 23 is a generalization of the Perron- Frobenius Theorem 
in this infinite-dimensional framework. 



Corollary 24. For any bounded and measurable function f , we have 

Ptf =l2(m) Y1 ^^^'^^/^ 



Proof. Fix t > and let / be a bounded measurable function on 
that Pff belongs to lJ{fi). On the one hand, we have 



(44) 

Let us first prove 



(P,/(x))Xx) < 



< oo. 



On the other hand, by Proposition 21, we have, for all x G M.*^, 
Ptfix) < 
< 

But the function 



^Qix)+lCt 



P 2 



lt<To(.)e-5«(-*) 



27rt 



dy. 



?/ I— J- e 2 



\Q(y) 



1 



-e 4 



is integrable on ]0, -|- oo[. Since e 2^(2/ < 1, we deduce that there exists a constant 
Kt> Q independent of x and / such that 



Ptf{x)<Kt\\f\\^e-^ 



and thus 



2 

00 ■ 



{Ptf{x)fd^^{x)<K^,\ 
Finally {PtfY integrable with respect to /i, so that Ptf G L^(/i). 
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Now we deduce from Theorem 1231 that 

Ptf =l2(m) ^{Ptf,r]i}^ r]i 



(45) 



ieN* 



If / belongs to IL^(/i), then the symmetry of Pt imphes that 

Since rji G IL^(yu), we deduce from the Dominated Convergence Theorem that the equahty 
{Ptf,r]i)^ = e^'^'* {f\Vi)tJ. extends to all measurable bounded functions. This and the 
equality ( 45 ) allow us to conclude the proof of Corollary 24 □ 



5.2.3 Existence of the Yaglom limit 



By Corollary 24, we have for any bounded and measurable function /, 



i>2 

< g-2(t-l)(A2-A 



Using Cauchy-Schwartz inequality, we deduce that, for any function h G IL^(/i), 



(46) 



By Theorem 22, S^Pi has the density r(l,x,.) G IL>^(/i) with respect to /i, so that, taking 
h = r(l,x,-), 



-.Alt 



^*P,+i/(x) - {v„f){v„r{l,x,-)),\ < e-2(*-i)(^-^^)||Pi/|k2(^)||r(l,x,-)||L2(^). 
By definition of (?7i, r(l,x,-))^ = e~^^rii{x). Thus we have 

e^^*P,+i/(x) > {m,f),e-'^r],ix) 

t— s>+oo 

and 

Finally, rii{x) being positive, for any x G M!j_, 

Ptfix) ^ ivij). 



«(/), 



where a is defined in (41 ). We conclude that a is the Yaglom limit for Z. We also deduce 



parts (2) and (3) of Theorem 20 
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5.2.4 Attractiveness of any initial distribution 

Let us first consider a compactly supported probability measure v on (0, + oo). By 



22 



y H-> J^, r{\^x^y)v[dx) is the density of vP^ with respect to [i. By 



Theorem 

Lemma 5.3]7 there exists a locally bounded function such that 

< 0(a;)?7i(y), Vx,?/ G (0, + oo) 



In particular, h : y ^ J^, r{l,x,y)i'{dx) belongs to L^. Then we deduce from (46) that 



To >t+i)= p^*^;y\ --^ «(/). 

We conclude that a attracts any compactly supported probability measure. 



Let us now prove that a attracts all initial distributions u supported in (0,oo). We want 
to show that, for any probability measure u on M^, for any Borel set A, we get 

lim F^{Xt e A\To > t) = a{A). (47) 



This is part (4) of Theorem 20 and it clearly implies the uniqueness of the QSD for X 
(and hence for Z). 

Proposition 25. For any a > 0, there exists ya > such that sup^.^^^^ E,x{e°''^y°-) < oo. 

Proof. Let us remark that e^^^^ e~'^^^^ dz dy < oo. Let a > 0, and pick Xa large 
enough so that e*^*^^-* ^-Qi^) dzdx < ^ . Let J be the nonnegative increasing func- 

" Xq^ <J X ACL 

tion defined on [xa,oo) by 



oo 



J{x) = / e«(^) / e-^^''Uzdy. 



Then we check that J" = 2g J' — 1, so that LJ = —1/2. Set now ya = 1 + Xa, and consider 
a large M > x. Ito's formula gives 

E^(^eaitAn,ATyj J(XmTmat,J) = J{x) + U {aJ{Xs) + LJ{X,)) ds\ . 

But LJ = -1/2, and J(X,) < J(oo) < l/(2a) for any s < Ty^, so that 

E^^,a(MTMAT,J J^X,,T,^^TyJ) < J{x). 

For X > ya, one gets l/(2a) > J{x) > J{ya) > 0. It follows that E^(e"(*^^'>^^^^«)) < 
l/{2aJ{ya)). Letting M -)■ oo then t oo , we deduce Ei:(e"'^"«) < l/{2aJ{ya)), by the 



monotone convergence theorem. So Proposition 25 is proved. □ 



Proving that a attracts all initial distributions requires the following estimates near and 
oo. 
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Lemma 26. For h G ^^{fJ') strictly positive on (0,oo) we have 

f^h{x)F,{To>t)fi{dx) 



lim lim sup ° = 0, 

lim limsup = 0. 

t^oo J^°° h{x)F^{To>t)fi{dx) 



(48) 
(49) 



Proof. We start with (48). Using Harnack's inequality (see [55| Theorem 1.1]), we have 
for £ < 1 and large t 

/i(x)P,(To > t)iJ,idx) ^ Pi (To > t) /i(^)/i(t/;2) 



then 



h{x)F,{To > t)fi{dx) C h{x)fi{dx)F^{To >t-l 



J^h{x)F,{To>t)fxidx) ^ F^iTo > t) h{z)fiidz) 

lim sup p , , , — ^ — , , , < lim sup 577; 

t-,00 Jh{x)F,{To>t)i2{dx) t^oo C/f '/i(x)/i(rfa:)Pi(To >t-l) 

e-^^ j;h{z)^i{dz) 

C J^^"^ h{x)fi{dx) 

and the first assertion of the lemma is proved. 

For the second limit, we set Aq := sup E^(e ^ ) < 00, where yx-^ is taken from Propo- 
sition 25 , Then for large M > yx^ , we have 

P,(To >t)= [ P,„(To > u)F,{T,, ed{t-u))+ P,(T,, > t). 
Jo 

Using lim e^^'^F^g{To > u) = ?7i(xo)(?7i,l)^, we obtain Bq := sup e^i"P^.o (Tq > u) < oc. 



u>0 



Then 



F..(To > t) < Bo [ e-^^"P,(T,„ e d{t - u)) + P,(T,„ > t) 
Jo 

< Bo e^^i* E^(e^i^^o) + g^^^* E^(e^^^^o) < e"^i*Ao(5o + 1) 



and (49) follows immediately (since x > Xq > yx^ ^ T^g < Ty^^). 



□ 



Let V be any fixed probability distribution whose support is contained in (0,oo). We must 



prove (47). We begin by claiming that v can be assumed to have a strictly positive density 



with respect to /i. Indeed, let 



+00 



^{y) = I r{l,x,y)iy{dx). 
Jo 



00 P + CX3 



JO 



+00 p+oo 



Using Tonelli's theorem we have 

r{l,x,y)u{dx) fi{dy) = I I r{l,x,y) fi{dy) u{dx) 

Jo Jo 

"+00 



r+oo 

/ p,(ro > i)u{dx) < 1, 
^0 



45 



which imphes that J r{l,x,y)i'{dx) is finite dy—a.s.. Finally, define h = i/ J idfi. Notice 
that for dp = hdfi 

F,{Xt+i G ■ I To > t + 1) = Fp{Xt G ■ I To > t), 

showing the claim. 

Consider M > e > and any Borel set A included in (0,oo). Then 

jF^iXt eA,To> t)h{x)p{dx) CF,{Xt e ATo > t)h{x)ii{dx) 



f P,(To > t)h{x) p{dx) p^(To > t)h{x) fi{dx) 

is bounded by the sum of the following two terms 



n 



12 



J F^iXt eA,To> t)h{x) fiidx) F,{Xt eA,To> t)h{x) p{d: 



x] 



jF^{To>t)h{x)fi{dx) 



jF^{To>t)h{x)fx{dx) 



F,{Xt e A, To > t)h{x) pidx) F,{Xt eA,To> t)h{x) p{d 



M . 



l''F,{To>t)h{x)fi{dx) 



/P^(To > t)h{x)fi{dx) 
We have the bound 

P,(To > t)h{x) fi{dx) + P,(To > t)h{x) fi{dx) 



Jl V /2 < 



/P^(To > t)h{x)n{dx) 
Thus, from Lemma [26] we get 

jF,{XteA,To>t)h{x)fi{dx) 



lim lim sup 

£4,0, MtoO 



jF,{To>t)h{x)pidx) 



l''F,{XteA,n>t)h{x)fi{da 



jfF,{To>t)h{x)fi{dx) 



0. 



On the other hand we have 



. //'P.(Xj G A,To > t)h{x)fi{dx) _ J^r],{z)fi{dz) 



lim ^ 

t—^oo 



aiA), 



f^^'F,{To > t)h{x)fi{dx) J^^Viiz)Kdz) 
since a attracts any compactly supported probability measures , and the result follows. 



Example 4. We develop now a numerical illustration of this logistic Feller diffusion case. 
As for the logistic birth and death process (see Example 3, Section |4]), the value of the 
charge capacity ^ will remain equal to the fixed value 9 across the whole example. 

We begin by showing in Figure [7] a random path of a logistic Feller diffusion process with 
initial size Zq = 1 and with parameters r = 9 and c = 1 (an Euler method is used for the 
numerical simulation of the random path). We observe that the process quickly attains 
the value of the charge capacity and remains around it for a long time. 

We compare now the Yaglom limits of two different logistic Feller diffusion processes 
whose charge capacity is equal to 9 (see Figure Isj): 



46 




Figure 7: Example 4. A random path for logistic Feller diffusion process with initial size Zq ~ 1 and 
parameters r = 9 and c = 1 



(a) 



(b) 



Figure 8: Example 4. The Yaglom limits of two logistic Feller diffusion processes with the same charge 
capacity: (a) r = 9 and c — 1; (b) r = 3 and c — 1/3. 
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(a) Z^"'\ whose parameters are r = 9 and c = 1, 

(b) Z^''\ whose parameters are r = 3 and c = 1/3. 

As for the logistic BD-processes, we observe that the two Yaglom hmits are centered 
around the charge capacity. But as a consequence of the relatively weak noise around the 
charge capacity, the Yaglom limit has clearly a smaller variation around this value in the 
logistic Feller diffusion case than in the logistic BD process case. We also observe that the 
smaller are the parameters, the flatter is the Yaglom limit and with a similar explanation 
as in the logistic BD-process case. 

We observe now the distance between the conditional distributions of Z^"^ and Z^^^ and 
their respective Yaglom limits, for different initial states, namely Zq = 1, Zq = 10 and 
Zq = 100. The results, computed with the help of the approximation method studied in 
Section [6| are represented on figure § For both Z^^^ and Z^^\ the speed of convergence 
to the Yaglom limit is the highest for Zq = 10, which is quite intuitive since the value of 
the charge capacity is 9. We also observe that it is higher for the processes starting from 
100 than for the processes starting from 1. In particular, this behavior is different than 
in the logistic birth and death process case. 
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5.3 The Q-process 

Let us now describe the law of the trajectories conditioned to never attain 0. 

Theorem 27. fTSf Let us fix a time s and consider B a measurable subset o/C([0,s],]R+). 
Then for any x G M.*^, 

hm P^(X G B\t < To) = Q^{X e B), 

t—^ca 

where is the law of a continuous process with transition probabilities given by q{s,x,y)dy, 
with 

q{s,x,y) = e^'' r{s,x,y) e-^^y\ 

r]i{x) 

Proof. Since r{s,x,y) e~'^^y^dy is the law of Xg started from x before extinction, we have 
to prove that 



7]i{X} 

For t > s, 

F,{XeB-To>t) ^ F^jX eB-To> s-ExATo>t-s)) 
Px(ro > t) P,(To > t) 

and we have proved that 

li^ ^yjTo >t-s) ^ rjM 
t^^ P,(To>t) r/i(x)' 

Then, 

t^oo P^(To > t) rii[x) 

Corollary 28. For any Borel set A C (0,oo) and any x, 

lim Q^{Xs e A) = r]l{y)fi{dy) =< r]i, 1 >^ / T]i{y)a{dy). 
Proof. Since 1a ^71 ^ IL'^(yu), thus 

r]i{x) Q^{Xs eA)= lA{y) Viiy) e^'' r{s,x,y) ^i{dy) 



□ 



converges to rji{x) f^ril{y)fi{dy) as s — )■ +oo, since e'^^'^r(s,x,.) converges to rii{x) rii{.) in 
lJ(dfi). □ 

Remark 8. The stationary measure of the Q-process is absolutely continuous with respect 
to a, with Radon- Nikodym derivative < ?7i, 1 >^ rji. 
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5.4 The case of a multi-type population 



Until now, we have considered a population where all individuals have the same ecological 
parameters. This biological assumption corresponds to the case where individuals have 
the same type. In this section, we generalize the previous study to a population composed 
of k different types. The population size process describing the dynamics of each subpop- 
ulation is given by a A;-dimensional stochastic Lotka-Volterra process Z = {Z}, ■ • • , Z^)t>o 
(SLVP), which describes the size of a A;- types density dependent population. This model 
generalizes to k types the 2-types density dependent model introduced by Cattiaux and 
Meleard 



More precisely, we consider for i,j G {1, ■ ■ ■ ,k} the coefficients 

7i > , > ; Cij> 0, Vi,j G {1, ■ ■ ■ ,k}. 
The process Z takes its values in and is solution of the stochastic differential system 



dZi = V^tdBl + (r.Zl - J2 ^vZlZi) dt, (50) 
where (i?*)j=i^... are independent standard Brownian motions independent of the initial 



data Zq. The system (50) can be obtained as (31) as approximation of renormalized 
fc-types birth and death processes in case of large population and small life lengths and 
reproduction times. The coefficients rj are the asymptotic growth rates of i-type's popula- 
tions. The positive coefficients 7j can be interpreted as demographic parameters describing 
the ecological timescale. The coefficient Cij, for i,j = 1, ■ ■ ■ ,k, represents the pressure felt 
by an individual holding type i from an individual with type j. Intra-specific competition 
is modeled by the rates cu, while inter-specific competition is described by the coefficients 
Cij > 0,i ^ j. If Cij = for all i ^ j, the stochastic /c- dimensional process reduces to k 
independent Feller logistic diffusion processes. Extinction of the population is modeled by 
the absorbing state (0, ■ ■ ■ ,0) and the extinction of the subpopulation of type i is modeled 
by the absorbing set 

Hi = {RlY^^ X {0} X {Rl)^-\ 



We denote by D the open subset of M'^ defined hj D = (M^)*^ and by dD its boundary. 
We denote by Tq the first hitting time of (0, ■ ■ ■ ,0), by T4 the first hitting time of some 
subset A and thus by Tqd the exit time of D. Of course, some of these stopping times 
are comparable. For example if the initial condition belongs to D, 

T9d<Th,<To, V2 = l,---,fc. (51) 

On the other hand, T/^. and T^j are not directly comparable for i ^ j. 

Let us prove the existence of the SLVP . 

Proposition 29. The process {Zt)t is well defined on ]R_|_. In addition, for all x G (1^+)*^, 

p,(ro < +00) = 1 

and there exists A > such that 

sup E^(e^^°) < +00. 

a:e(K+)'= 
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Proof. The existence of the SLVP is shown by a comparison argument (cf. Ikeda-Watanabe 
[36j Chapter 6 Thm 1.1). Indeed, the coordinates {Z^t can be upper-bounded by the in- 
dependent solutions of logistic Feller equations 

dV; = ^/^^dBl + {r,Yl - cuiYlf) dt, (52) 

for which we have obtained in the previous section that extinction occurs a.s. in finite 
time and that the extinction time has some finite exponential moments. The almost sure 
finiteness of each T^-., hence of Tqd and Tq, thus follows. □ 



As in the previous sections we are interested in the quasi-stationary distributions for the 
process (50). We firstly reduce the problem by a change of variable. Let us define Xt = 



. Xn with XI 



We obtain via Ito's formula and for any i G {1, ■ ■ ■ ,k}, 



«; ^ .Bi.llf -^S^2iM--L), (53, 



2X? 



In the following, we will focus on the symmetric case where X is a Kolmogorov diffusion, 
that is a Brownian motion with a drift in gradient form as 

dXt = dBt - VV{Xt)dt. (54) 



Let us state a necessary and sufficient condition to write the drift of (X) as in (54). The 
proof is obtained by computation and requires the equality of the second order cross- 
derivatives of V. 

Proposition 30. // the following balance conditions on the ecological parameters are 
satisfied, 

Cijlj = Cjai, Wi,j, (55) 
then the process X is a Kolmogorov process with potential V given by 

v{x\ ■ ■ ■ ,x^) = I + '-""^ - ^) + E^^.^^- (^^)' (^'■)'- 

i=l ^ ^ ij^j 



We will establish an existence and uniqueness result for the QSD of the process (Xt). 
The re-statement of the results for the initial stochastic Lotka-Volterra process follows 
immediately, since the hitting time of (0, ■ ■ ■ ,0) and Hi and the exit time of D are the 
same for both processes (X) and (Z). 

By generalizing to A;-types populations the results proved in [T3] for two-types populations 
(an easy consequence of Girsanov's theorem), we get 

Proposition 31. For all x & D, for all i ^ j , 
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Let us now state the first theorem, which is concerned by conditioning on the co-existence 
of the k types. 



Theorem 32. Under the balance conditions (55), there exists a unique quasi-stationary 
distribution v for the process {X) and the absorbing set dD, which is the quasi-limiting 
distribution starting from any initial distribution: for any ^ on D and any A (Z D, 



lim F^{Xt G AITqd >t) = u{A). 

t—^+oo 



Furthermore, there exist A > and a positive function rj such that 

lim e^'F^iXt e A; Tqd > t) = r]{x) u{A). 

t— s>+oo 

Proof. The proof of the existence of a quasi-stationary distribution results from the spec- 
tral theory for the semi-group of the killed process (Pj) (related to X) established in 
Cattiaux-Meleard [H], Appendices A, B, C. Define the reference measure on (M+)^ by 

fi{dxi, ■ ■ ■ ,dxk) = e"^^*-^^ dxi ■ ■ ■ dxk- 

As in Subsection 5.2.2, one builds a self-adjoint operator on L^(/i) which coincides with 
Pt for bounded functions belonging to L^(/i). Its generator L is self-adjoint on L^(/i) and 

Lg = ]^/\g-V -Vg, \/gEC^iD). 

We check that the assumptions required in [H] Theorem A. 4 are satisfied and therefore, 
the operator — L is proved to have a purely discrete spectrum of non-negative eigenvalues 
and the smallest one A is positive. The corresponding eigenfunction t] is proved to be in 
L^(;u) and the probability measure u = j^^^^ is the Yaglom limit. 

let us emphasize that the uniqueness of the quasi- stationary distribution results by [14] 
Proposition B.12 from the ultracontractivity of the semi-group Pt (ultracontractivity 
means that Pt maps continuously L^(yu) in L°^(/i) for any t > 0). The proof of the 
latter is easily generalized from the two-types case ([H] Proposition B.14) to the A;-types 
case. 

□ 



Theorem 32 shows that in some cases, a stabilization of the process with co-existence of 
the k types will occur before one of these types disappears. Let us now come back to our 
initial question: the long-time behavior of the process conditioned on non-extinction. For 
each i = 1, . . . ,k, we denote by Aj the smallest eigenvalue related to the purely discrete 
spectrum of the generator for the i-axis diffusion defined by the stochastic differential 



equation (52). 



Theorem 33. Under the balance conditions (55), there exists a Yaglom limit m for the 



process (X) conditioned on non extinction: for any x ^ 0, for any A C D, 

lim P^(Xt G A\To >t) = m{A). 

The support of this measure is included in the k axes. 

Furthermore, if there exist ii,...,ii G {1, ■ ■ ■ ,k} such that A*^ = • • • = A*' < minj^jj^...^j; A*, 
then this QSD is concentrated on the axes of coordinates ii,...,ii. 
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Proof. Recall that the existence of a Yaglom limit has been proved in the case k 



(Section 5.2). In what follows, we prove by induction the existence of a Yaglom limit for 



any A;-type system (53) 



The induction assumption (Afc_i) is as follows: we assume that, for any {k — l)-type 



Kolmogorov process X^'^ satisfying (53) with (55), there exist a constant A > 0, a 
uniformly bounded function rj > and a probability measure u on (M+)'=-^ such that, 
for any x G (]R+)'^^^ \ {0} and any bounded measurable function / on (]R+)*''^^ such that 
/(O) = 0, we have 



hm e^*E,.(/(xf-^)))=r/(x)K/); 



t—^co 



sup 

t>o, xeCM*.)*^-! 



|e^*E,.(/(xf-^)))|<+oo. 



(56) 



As mentioned above. Assumption (Ai) is already proved. Let us assume that {A^^i) is 
true and show that (Ak) follows. 

Let X*^*"'^ be a A;-type Kolmogorov process satisfying (53) with (55). Once hitting the 



boundary dD = U^^^ifj, the process will no more leave it. Hence, for t > Tqd, the 

process will stay on the union of hyperplanes Hi. Moreover, Tqd = inf ^ Tji^^. Fix 

i e {1, ■ ■ ■ ,k} and assume that the process leaves D through Hi. The dynamics on is 
given by the process iui'^'^)j^i defined in (M+)^-i by: 



da 



dBl + 



2U, 



Remark that by Proposition 31, the process really leaves dD by the interior of H^. Each 



system {U^''^'^)j^i is a (A; — l)-type kolmogorov process (53) with balance conditions. Hence, 
by our induction assumption (Ak-i), there exist for each i G {1, ■ ■ ■ ,k} a positive constant 
Vi, a positive function rji and a probability measure z/j on Hi such that (56) holds for 
(f/«'^W„zG{l,---,fc}. 

Let us define 

Vmin = inf Vi. 

ie{i,-,k} 

For any bounded measurable function / on (M+)'^' such that /(O) = and for all t > 0, 
we have 



i=l 



t )^Tgn=TH,<t ) 5 



(57) 

where we used the fact that X*^^^ reaches dD by hitting the interior of one and only one Hi. 



By Theorem 32, there exist a positive constant A', a positive function rj' and a probability 
i*,)'' such that 



measure u on 



hm e^'%(/(xf))lr,,>,) =V(x)z.'(/). 
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Moreover, a similar coupling argument as in [H] yields A' > fmm- We deduce that 

lim e^— % (/(Xf )lT,,>i) =0. 
For each i G {I,-- - ,k}, we have by the Markov property 



By the induction assumption (Ak-i), e''™"*^*~"^'3°'*/(f//l''j,^^) is uniformly bounded. More- 
over the inequality < f , < A' and Proposition [s] ensure that E^. (^e^™'"^^-°) < +oo. 
Using the convergence property of the induction assumption {A^-i) and the dominated 
convergence theorem, we deduce that 

lim E. (e-^Vlxf ^)1t,,=t..<.) = | (e"""'^^^lT..=T.^r].(X^,j) ^.(/), if v. = v^.^ 
^ ' ^ 1 0, otherwise. 



We have then 

k 

hm e^-*E,(/(xi'=))) = 5^U,=,_E,. (e^--^^-lT,,=T„M^Ten)) Mf), 

i=l 



which gives us the first part of the induction assumption (Ak). 

In order to prove the second part of (Ak), let us introduce the SLVP F^^-* with coefficients 
{c'ij) defined by 

c'kk = Ckk, Aj = Cij and 4^ = c^^ = 0, Vz,j = 1, ■ ■ ■ ,k - 1. 

By the same coupling argument as above, the return time to dD for X^''^ is stochastically 
dominated by the return time to dD for Y^''\ i.e. P^(xf ^ E D) < P^(f/''^ G D) for all 
t > 0. 

Since the k — 1 first components of Y^'^^ are independent of the last one and since 
we have 

p,(r/^'^ G < p,.((f/''^'\ ■ ■ ■ g x p,(r/^)''^ g m;). 

On the one hand, the dynamic of ■ ■ ■ ^yC^).^-!) is the same as U^''\ so that, by the 

second part of the induction assumption {A^-i) and by the definition of Vmin, 

sup e"— *P,((rf'^'\ ■ ■ ■ ^''^-i) G (M;)'-^) < +CX). 
t>o,a;ei:> 

On the other hand, y(^)''^ is a one dimensional SLVP, thus we deduce from (Ai) that there 
exists a positive constant Ai such that 

sup e^^*P^(r/*^^''' G M;) < +00. 
t>o,xe-D 
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As a consequence, we have 

sup e(''™'"+^^)*P^(xf ^ eD)< sup e^'"rmn+Xi)t^^(^^lk) ^ < 

t>0,x£D t>0,x£D 

and we deduce that 

supE^(e''™'"'^^^) < +00. 
For any bounded measurable function /, this immediately leads us to 

sup E,(e^— *lt<T,,/(xf )) < +00. 

t>0,xGE 

Moreover, by Equality (58) and the second part of (A^^i), we deduce that, for each 
ie {!,■■■ ,k}, 

sup (e^-V(xf ^)lTa.=T,^<*) < 
t>o,xeE ^ ' 



-00. 



By Equality (57), the second part of the induction assumption (Afc) is thus proved. 

By induction on > 1, we conclude that Assumption (A^) is true for any /c > 1, thus 
Theorem [33] follows. □ 

Example 5. Let us numerically study a 3-type system and observe its long-time behavior. 
The 3-tuple process (Z^, Z^, Z^) evolves as 

dZ\ = ^f^ldBl + [r^Zl - c^^{Zlf - c^^Z\Zl - c^^Z\Zl) dt, 
dZ^ = y^dB^ + {r,Z^ - c^iZ/Z^ - c22iZff - c^sZ^Zf) dt, 



with 



and 



dZ^ = V^dBf + {r,Z^ - c,^ZlZ^ - c.^Z^Z^ - c,,{Zlf) dt, 



7i = l, Cii = 10, V^G {1,2,3} and c,,- = 0.5, V 2 7^ j G {1, 2, 3}, 



ri = 1.5, r2 = 1, rs = 0.5 ; Z] = Zl = Z^ = 1. 



We describe the dynamics of F(^i^i^i-j{{Z^ ,Z^ ,Zf) G ■ | Tq > t). As explained above, the 
process conditioned on non-extinction initially behaves as a 3-type population. Then a 
type goes extinct, then a second one and finally it only remains one type in the population. 
In order to represent graphically these transitions, we compute numerically the dynamics 
of the probabilities of coexistence and existence of the different types as functions of time. 
In Figure 10 , we represent 

(a) the probability of coexistence of the three types P(i 1 1) {Z^ > 0, Z"^ > 0, Z| > | Tq > 

t); 

(b) the probability P(i > 0, > 0, Z^ = \ Tq > t) of coexistence of exactly two 
types i 7^ j, for each combination of types {i,j,k) = (1,2,3), {i,j,k) = (2,3,1) and 
ii,j,k) = (1,3,2); 
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^(1,1,1) (2j'>o,z,^>o,z;'>o|ro>t) 



(a) 



,1,1) {Zl>O.Z^>O.Z^^O\TB>t) 




To > t] 



0\T„>f] 



(b) 



P(i,i,i) {Zl >(i,Zf = 0,Zf =0\To>t 




Figure 10: Dynamics of the probabilities of co-existence and existence of the different types for a 3-type 
stochastic Lotka-Volterra system. The horizontal axis is the time axis 
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Figure 11: First marginal of the Yaglom limit of a 3-type stochastic Lotka-Volterra system. The two 
other marginals are equal to the null measure. 



(c) the probabihty P(i^i^i)(Zj > 0, = 0, = | Tq > t) of existence of one and only 
one type i, for each type i = 1,2 and 3. 

As expected, the 3-type mode disappears quickly and the 2-type modes are transient. 
We also observe that the probability P(i,i,i)(Z/ > 0, = 0, Z^^ = | Tq > t) con- 
verges to 1 when t increases, meaning that the last state of the population before ex- 
tinction is monotype with type 1. It turns out that the support of the conditional law 
F(^i^ii-^{{Z},Z^,Z^) G ■ I To > t) becomes more and more concentrated on x {0} x {0} 
in the long time. The Yaglom limit is thus equal to Ui (g) 5(o,o)) where Ui is the Yaglom 
limit of the process 



dZl' = VZi'dBl + {t,Z[' - Cn(Zf )2) dt, 
absorbed at and is represented in Figure 11 



6 Simulation: the Fleming- Viot system 

As seen in the previous sections, the spectral theory is a powerful tool to prove existence 
and eventually uniqueness of a QSD for a given process Z. It is based on the equivalence 
property of Proposition [4], stating that a probability measure a on E* is a QSD for the 
killed process Z if and only if 

aL = —9{a)a, (59) 

where L denotes the infinitesimal generator of Z and 6 (a) a positive constant. In some 
cases, such as in the finite state space case, one can easily compute numerically the whole 
set of eigenvalues and eigenvectors of L as seen in Example 1 and Example 2. For these 
numerical illustrations, we used the software SCILAB and its function spec. We also 
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refer to [62] for a detailed description of some algorithms available in MATLAB for the 
computation of eigenfunctions and eigenvalues in large (but finite) state space cases. 

In other cases, such as the logistic birth and death process of Section |4] and the logistic 



Feller diffusion of Section |5| solving numerically Equation (59) is too hard and we use a 
different approach. This approach consists in approximating the QSD and the conditioned 
distribution Fz{Zt E .\t < Tq) by the empirical distribution of a simulable interacting 
particle system. This Fleming- Viot type system, built for any number of particles N > 
2, has been introduced by Burdzy, Holyst and March jlO] and explored in [11] and in 
Grigorescu-Kang [32j for rf-dimensional killed Brownian motions. It has also been studied 
in Villemonais [65j for multi-dimensional diffusion processes with unbounded drifts and a 
general result is available in [SB]. Similar systems have also been considered by Ferrari- 
Marie [23] for continuous Markov chains in a countable state space. In this section, we 
explain the approximation method based on the Fleming- Viot type interacting particle 
systems. 

Let Z he a killed Markov process which evolves in the state space E. Fix N > 2 and let 
Zq e £" be its initial value. The interacting particle system with particles {Z^, ■ ■ ■ ,Z^) 
starts from (Zq, ■ ■ ■ ,Zq) and belongs to [E*)^ . The particles evolve independently from 
this initial position according to the law of the killed Markov process Z, until one of them 
hits the state 0. At that time ri, the killed particle jumps to the position at ri of one of 
the A^ — 1 remaining particles, chosen uniformly among them. Then the particles evolve 
independently according to the law of Z until one of them attains (time r2), and so on. 
The sequence of jumps is denoted by (r„)„ and we set 

Too = lim r„. 

n— >oo 

This procedure defines the (ii^*) ^-valued process (Z^, ■ ■ ■ ,2'^) for all time tG [0,roo[. 
Figure [T2] shows an illustration of such a system with two particles evolving between their 
jumps as Markov processes absorbed in and 1. 

If Too = +00 almost surely, then the Fleming- Viot particle system will be well defined at 
all time t > 0. The condition Too = +oo is clearly fulfilled for continuous time Markov 
chains with bounded jump rates. In the diffusion process case, criteria have been provided 
in [g, [33], [65] and [66]. 

In that case, denote by the empirical distribution of (Z^, ■ ■ ■ ,Z^) at time t: 

1 ^ 

i=l 

The following result is obtained in [66] by martingale method. 

Theorem 34. Assume that for all N > 2, {Z^, • • • , Z^) is well defined at any time t >0. 
Then, for any time t > 0, the sequence of empirical distributions {^f) converges in law 
to the conditioned distribution P^^ [Zt E - {t < Tq), when N goes to infinity. 



If moreover {Z , ■ ■ ■ , Z ) is ergodic, we denote by M its stationary distribution and by 

N 



its empirical stationary distribution, which is defined by = -^"^f^iSz,, where 
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"''1 T2 "^4 

Figure 12: Fleming- Viot type system with two particles absorbed in and 1. 



{zi, ■ ■ ■ ,zn) & E* is a random vector distributed with respect to M^. In particular, fi^ 
converges in law to X'^ when t — )■ oo. We refer to [65] for the proof of the following 
theorem. 

Theorem 35. Assume that Z has a QLD a which attracts all initial distributions: for 
any probability measure fi on E* , 

lim ¥^{Zt e ■\t<To) = a. 

i— s>+oo 

Assume moreover that {Z^, ■ ■ ■ ,Z^) is ergodic and that the family of laws of (A'^)jv>2 is 
uniformly tight. Then the sequence of random probability measures {X'^) converges weakly 
to a. 

If ii^ is a bounded subset of M'^, d > 1, and if Z is a drifted Brownian motion with 
bounded drift which is killed at the boundaries of E, then the assumptions of Theorems 



34 and 35 are fulfilled (see |65)). The proofs of Theorems 34 and 35 are based on a 
coupling argument. More general (but longer) proofs can also be found in ||33| or |66| . 
In particular, these results provide us a numerical approximation method of the Yaglom 
limit for such processes. 



Let us now consider the Kolmogorov diffusion process X defined in (39). In that case, 
the existence of the Fleming- Viot particle system remains an open problem because of 
the unboundedness of the drift coefficient. In order to avoid this difficulty, we introduce 
the law of the diffusion process with bounded coefficients defined by 

dX^ = dBt - q{Xl)dt ■ Xo G {e, 1/e), (60) 
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killed when it hits e or -. One can easily show that at any time t > 0, the conditioned 
distribution of converges to the one of X: 

r (Xf e-\t<T,A Ty,) — -> P (Xt G -It < To) . 

The existence of the Yaglom limit denoted by and the uniqueness of the QSD for 
are obtained from Pinsky [49j . The following approximation result is proved in [65J using 
a compactness- uniqueness argument. 

Proposition 36. The sequence (a^)^ weakly converges to the Yaglom limit a of X as e 
tends to 0. 



For all X > 2, we denote by (X^'-"^,-- - ,X^'^) the interacting particle system built as 
above, with the law P^. Since the diffusion process X'^ is a drifted Brownian motion with 
bounded drift evolving in the bounded interval ]e,l/e[, the interacting particle system 



(X^'^, ■ ■ ■ ,X^'^) fulfills the assumptions of Theorems 34 and 35, Denoting by /x*^'^ the 



empirical distribution of the simulable particle system (X^'^, ■ ■ ■ ,X^'^), we get 



and 



lim lim /i^'^ = Pxn {Xt € -It < Tq) , Vt > 

e->0 N^oo 



lim lim lim /i^'^ = lima^ = a. 

e— !>0 N-^oc t— s>oo e->0 



Then, choosing e small enough and X big enough, we get a numerical approximation 
method for the conditioned distribution and the Yaglom limit of X. 

Example 6. Let us now develop this simulation method in the case of the Wright-Fisher 
diffusion conditioned to be absorbed at 0, which evolves in [0,1[ and is defined by 



dZt = ^Zt{\ - Zt)dBt - Ztdt, Zo = z G]0,1[. 

This is a model for a bi-type population in which the second type cannot disappear. In 
that model, Zf is the proportion of the first type in the population at time t > and 
1 — Zf the proportion of the other one. The existence of a Yaglom limit for this process 
has been proved by Huillet in [35], which also proved that it has the density 2 — 2x with 
respect to the Lebesgue measure. 

Using the approximation method described above with e = 0.001 and X = 10000, we 



obtain numerically the density of the Yaglom limit for Z represented in Figure [13} which 
is very close to the function x t— )■ 2 — 2x and shows the efficiency of the method. 
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